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ABSTRACT 


A  method  for  the  computation  of  confidence  intervals  for 
circular  error  probability  (CEP)  based  on  first  order  variance 
estimates  was  introduced  in  1966.  It  was  later  found  that 
under  certain  conditions  the  resulting  confidence  intervals 
for  CEP  were  smaller  than  expected.  As  a  result  a  second 
order  variance  estimate  method  was  developed,  at  the  Johns 
Hopkins  University  Applied  Physics  Laboratory,  which  greatly 
improved  the  accuracy  of  the  confidence  intervals  for  CEP. 
The  purpose  of  this  thesis  is  to  develop  and  test  procedures 
for  the  3-dimensional  case  to  obtain  a  second  order  estimate 
for  variance  of  spherical  error  probability  (SEP) . 


iii 


Acoeaslon  Pop 

X 

NTIS  cnAJtl 

DTi": 

□ 

I’l  ■  t.  ’nc"  tri'TiJd 

□ 

8y.-_ 

I  Distribution/ 

Availability  Co<!<'s 
lAvall  an<3/or 
Spocial 


THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed 
in  this  research  may  not  have  been  exercised  for  all  cases  of 
interest.  While  every  effort  has  been  made,  within  the  time 
available,  to  ensure  that  the  programs  are  free  of 
computational  and  logic  errors,  they  cannot  be  considered 
validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I .  INTRODUCTION 


A.  BACKGROUND 

In  1966  W.R.  Blischke  and  A.H.  Halpin  [Ref.  1]  first 
described  methods  for  the  computation  of  confidence  intervals 
for  circular  error  probability  (CEP)  based  on  first  order 
variance  estimates.  When  these  procedures  were  implemented  at 
the  Johns  Hopkins  University  Applied  Physics  Laboratory  (JHU- 
APL)  it  was  found  that  under  certain  conditions  the  resulting 
confidence  intervals  for  CEP  were  smaller  than  expected.  As 
a  result  R.C.  Ferguson,  K.V.  Kitzman  and  P.B.  Jackson  [Ref.  2] 
have  extended  the  Blischke~Halpin  procedures  tc  create  a 
second  order  variance  estimate.  Testing  conducted  by  Kitzman 
[Ref.  3]  has  confirmed  that  this  method  produces  a  more 
accurate  estimate  for  the  variance  of  CEP. 

Our  goal  will  be  to  extend  this  procedure  to  the  3- 
dimensional  case  to  obtain  a  second  order  estimate  for 
variance  of  spherical  error  probability  (SEP) . 

B.  METHODOLOGY 

We  begin  by  assuming  that  missile  detonations  are 
distributed  as  trivariate  Gaussian.  In  the  past,  analysis  of 
ballistic  missile  test  firing  data  has  failed  to  disprove  this 
assumption.  If  the  mean  and  variance  are  given  by: 


1 


=  \^2  . 


Then  the  probability  of  a  detonation  occurring  within  a 
distance  R  of  the  origin  is  given  by: 

=  fff  — A— exp{--|(;c-n)^-\jr-p)ldxr,dy,dx, 

•'V  (^/2Wfq  I  2  j 

where, 


g2  =  det[2], 

and  D  =  {(x^,  x^,  ^ 


By  the  implicit  function  theorem,  >  0  implies  that  the 

equation  P{R;\i,'L)  =  CONSTANT  defines  a  differentiable  function 
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such  that  S);  |i,  E)  =  CONSTANT.  If  the  CONSTANT  is 

set  to  ,  then  is  the  SEP  function.  Let  X^,  X.^, be 

a  random  sample  of  independent  and  identically  distributed 
random  vectors  with  distribution  .  Then  the  maximum 

likelihood  estimators  of  n  and  S  are  p  =  — 

£  =  —  ^  .  When  (1  and  £  replace  p  and  £  in  the 

equation  then  i?(il,£)  gives  SEP,  an  estimator  of 

SEP. 

The  second-order  approximation  for  the  variance  of  SEP  is 
given  by: 

*  -^(P:.„SEPl*'|f®P)fe.SEP) 

where  £“  =  (Oii ,  ,  0^2 , 0^3 ,  O33) .^uSEP  is  the  derivative  of 

SEP  with  respect  to  n,E"  viewed  as  a  column  vector,  D^^,.SEP  is 


3 


the  usual  second  derivative  matrix  (or  Hessian)  of  SEP  with 


respect  to 


is  the  variance-covariance 


matrix  ji,  ®  is  the  tensor  product  operation  and  the 

underline  notation  is  used  to  represent  the  column  vector 
formed  by  stringing  out  the  elements  of  the  matrix  by  rows. 


For  example, 


if 


A  = 


1  2 
3  4 


then  ^  =  [1  2  3  4]^. 


Note  that  jl  is 


distributed  so  P 

n 


It  can  also  be  seen  thatnU 


has  a  Wishart  distribution  with  parameter  S  ,  so  that 
“  (2(S*S)  .  The  complete  derivation  of  equation  (1)  ,  based 

on  Taylor  Series  expansion,  is  given  in  appendix  A. 
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since  (i  and  2  are  independent  [Ref.  4:p.  102]  it  follows 


Derivation  of  the  first  order  terms  in  equation  (2)  has 
been  accomplished  by  S.D.  Hill  at  JHU-APL  [Ref.  5],  In  order 
to  solve  for  the  second  order  terms  it  will  be  necessary  to 


find  expressions  for 


d^SEP 


d^SEP 

82“' 


d^SEP 

8^l82“ 


i.e.  it  is 


necessary  to  find  j,^SEP ,  the  second  derivative  matrix  of  SEP 

with  respect  to  its  arguments.  Chapter  II  will  describe  the 
derivation  of  the  secor.d  order  variance  estimate  and  Chapter 
III  will  discuss  implementation  and  testing  of  the  algorithm. 
Conclusions  will  be  presented  in  Chapter  IV. 
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II.  DERIVATION  OF  THE  SECOND  ORDER  ESTIMATE 


A.  FIRST  ORDER  TERMS 

We  begin  by  describing  the  derivation  of  the  first  order  terms 
given  in  equation  (2) .  If  we  define: 


where 


0^2 

0^2 

=  o2^ 

o21 

022 

023 

/ 

= 

022 

022 

023 

=  o22 

+(X2 2)2 0^2  +  2(^2 -|i2)(^3 -^A 3)0“ +(X3 -H3)2  0^3 


and 

g(p ,  0,  <|) ;  H ,  S“)  =  F(psin6cos4),  psin6sin(}),  pcos0;  p,  E“) . 


Then,  with  a  shift  to  spheiical  coordinates,  we  see  that: 
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We  wish  to  solve  equation  (3)  for  each  of  the  partials 


dSEP 

dX, 


where  A.  -  (^i/ ^2' ^*'3' °i2' ®i3' ®22' *^23' *^33)  •  From 


previous  section,  is  the  SEP  function,  so  we 


to  solve  for  the  partials 


dR 
dX^  * 


Taking  partials  of  equation  (3)  we  see  that. 


P'(R)  ^ 


31, 


and 


P'(R) 


sin8g(P,  0,  <J))d8d4) . 


We  can  solve  the  right  hand  side  of  equation  (4) 
i=l,2,3  by  using  the  relationships 

^9  --  ^9  3g  __  dg  3g  __  3g 

3nj^  dx^  d\L^  dx^  3^3  3^3 


the 


need 


(4) 


for 
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so  we  see  that 


If  we  apply  Green's  Theorem  to  this  equation  and  then  make 
a  shift  to  spherical  coordinates  we  get: 

-  j j  j dx^dx2dx2  =  /?^J^"J’'cos(t)sin^0g(i?,  0,  <t))d0d(j) 

n  1  0  0 


j j  j dx-^dx^dx^  =  /?^J^’'J"sin<|>sin^0g(i?,  0,  (t>)d0cf(J) 


-///^  dXjdXjdXj  =  i?^J''*J’'cos0sin0g(i?,  0, 4))d0d4> 

D  0  0  0 

=  f”sin26g(R,6,<i>)ddd<i> 

2  Jo  Jo 


To  solve  the  right  hand  side  of  equation  (4)  for  i=4 , .  .  .  ,  9 
we  first  note  that: 


dg 

ao,^  ' 

II 

dg 

ao,,  ' 

j*j 

8 


so  we  get 


^9 
do  • 


1  3 

2  8^^ 

_  3  Bg 
Bxj  d\i^ 


i 


which  is  used  to  evaluate  the  integral  of 


as  follows: 

81  _j 


+  o“(pcose-H3)]cos4)sin20gr(ie,  0,  <j))ded<t) 
'///  -§^^i^2^2  =  ^^f^'‘f''[<J^\psinQcos(\>-\x^)+a^\psin6sin^-\iA 

D  ^  ^  ^ 

-^a^3(pcos0-p3)]sin4)sin20g(i?,  0,  <j))d0d(}) 

+a“(pcos0-p3)]sin205t;?,  e,(t>)ded<p 

+a23(pcos0-p3)]sin<t>sin20g(/?,  0,  <j))d0d(t) 

~fff  "  ^/^'YoV>3^'^®^°3‘<>~Pi)*o“(Psin0sin<t>-p2) 

+o23(pcos0-p3)]sin20g(i?,  0 ,  (|>)c?0d<J> 


9 


If  we  define  the  functions: 

CC.  .  j(J?;  ji,  S“)  =  [*cos\jQ)cos\l^)g{R,Q,^)d6d^ 

2")  =  j'^’'J‘Jcos\jd)sin^(l<t>)g(R,d,<t>)ded(t> 

SC.  j(i?;  ji,  E“)  =  f^"f’‘sin^(j0)cos^(i4>)g(i?,0,4))d6d<|) 

SS^  •  =  r^’'r*sin'‘(j0)sin*(i<t>)g(-K»0»4>)c®c?4> 

'  ■  J  0  J  0 


Extending  the  work  of  Blischke  and  Halpin,  with  the  aid  of 
Green's  Theorem,  we  can  see  that: 


^^2 ,1,1,1 

^■^2,1,1,! 

■^^1 ,1,0,0 

^^1 ,1,0,0 

‘^'"1 ,2,0,0 
^  ,1,0,0 


T 


(5) 
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Applying  the  partials  given  in  equation  (6)  to  the 
partials  -JJJ'-^dx^dx2dXj  given  above  we  see  that: 
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dSEP 

do.~ 


dSEP 

^023 


dSEP 

8033 


where  Pv(;?;n,E‘^)  is  given  by: 


= 

,1,1,1 

*^2 

= 

^‘^1,3,1,! 

*^3 

= 

^•^1,1,1,! 

= 

^•^1 ,3,1,1 

*^5 

= 

■^^1 ,1,0,0 

>^6 

= 

■^^1 ,2,0,0 

= 

■^^1 ,3,0,0 

j‘^8 

= 

*^^2 ,1,1,1 

W, 

= 

*^^3 ,1,2,1 

*^10 

= 

^^2 ,1,1,1 

‘^11 

= 

■^■^3 ,1,1,2 

= 

■^■^3 ,1,2,1 

where  each  of  the  functions  is  evaluated  at  /?(n,E"),  n,  E“. 

This  completes  the  derivation  of  the  first  order  terms 
required  in  the  computation  of  the  variance  estimate. 
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B.  CONCEPT  FOR  DERIVATION  OF  THE  SECOND  ORDER  TERMS 

We  first  observe  that  equations  (5)  and  (6)  enable  one  to 

express  ^uSEP  as  a  composition  of  three  functions 

Tj  and  ,  each  of  which  is  a  fairly  simple  function  of  its 

arguments . 

Explicitly, 

=  r3°r2°T,{n,S“) 

where , 

E“). 
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since  ■^uSEPi^i , 'L'^)  =  we  may  apply  the  chain 

rule  to  obtain  the  second  derivative  as  a  matrix  product: 

Dl^.SEP{v^,i:-)  =  DT^iT^-T.iiX,  E“))-Z?r2(r>  ,  2“))-Dr>  ,  E“) .  (7) 


The  simplicity  of  the  functions  T^,  and  means  that 
the  derivatives  DT^,  DT^  and  DT^  are  easily  obtained  and  then 
the  rather  large  matrix  products  involved  in  DT^-DiyOT^  can  be 
left  to  the  computer  leading  to  a  rather  painless  computation 
of  D^  ^uSEP. 
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DT^  (the  total  derivative  of  the  function  )  is  the  10x9 
matrix  consisting  of  the  9x9  identity  matrix  (from  ) 

with  the  tenth  row  consisting  of  the  partials  of  SEP  with 
respect  to  mu  and  sigma. 

1  0  0  0  0  0  0  0  o' 

010000000 
001000000 
000100000 
000010000 
DT^  =0  0  0  0  0  1  0  0  0. 

000000100 
000000010 
000000001 
dSEP  dSEP  dSEP  dSEP  dSEP  dSEP  dSEP  dSEP  dSEP 
6^i3  ao,,  ao,2  ao^^  ao^j  ao33 


DT2  is  the  22x10  matrix  consisting  of  the  10x10  identity 
diu,  R) 

matrix  (from  — dt; - i — '),  with  rows  11  to  22  consisting  of  the 

,  2 

12x10  sub-matrix  of  partials  of  W  with  respect  to  ji,  2"  and 
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1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

dwj 

8W3 

8W3 

8W3 

8W3 

8W3 

8W3 

8W3 

a^A3 

8033 

8012 

8033 

8023 

8033 

dR 

dw^2 

8W32 

8W32 

8Wi2 

3*^12 

dw 

^'^12 

8W32 

d\i^ 

^^2 

8^3 

8011 

80i3 

^®22 

^023 

8033 

dR 

DT^  is  the  9x22  matrix  consisting  of  the  partials  of 
u^(i=l , ,9)  with  respect  to  n,  i?  and  W. 


■  8u, 

8u, 

du. 

du. 

-^(1x5) 

8S“ 

DT,  = 

8Uq 

8Uq 

8uo 

dU:, 

- ^(1x5) 

8S“ 

Thus,  in  order  to  obtain  the  second  degree  variance 
estimate  given  in  equation  (2)  we  must  solve  for  the  values  of 
all  elements  of  DT^,  DT2  and  DT^.  Once  these  values  are 

computed,  the  solution  is  found  by  the  matrix  multiplication 


16 


shown  in  equation  (7)  ,  yielding  the  9x9  matrix  of  second 
partials  of  SEP  with  respect  to  mu  and  sigma. 

The  values  of  the  DT^  elements  are  computed  from  equations 

(5)  and  (6)  .  The  determination  of  DT^  requires  the  evaluation 

dW  Sw 

of  ,  -  and  which  is  fully  described  in  appendix  B. 

Finally,  to  solve  for  the  elements  of  DT^  it  is  necessary  to 

solve  for  each  of  the  partials  given  above.  The  evaluation  of 
these  partials  is  given  in  appendix  C. 
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III.  IMPLEMENTATION  AND  TESTING 


A.  PROGRAMMING 

The  calculation  of  the  second  order  approximation  for 

variance  of  SEP  (when  |i  and  S  are  known)  has  been  implemented 

in  the  FORTRAN  program  SEPCMP  which  is  provided  in  appendix  D. 

The  following  steps  are  followed  to  obtain  the  approximation: 

Input  the  distribution  mean  (n)  and  variance- 
covariance  matrix  elements  (S“) . 

Compute  SEP  using  ^  and  S  . 

Compute  the  values  of  each  element  of  the  matrices 
DTI,  DT2  and  DT3  and  then  form  the  matrix  product 
comprising  the  second  derivative  of  SEP  with 
respect  to  ^  and  The  first  derivative  of  SEP 

with  respect  to  ^  and  is  then  given  by  row  10 
of  the  matrix  DTI. 

Calculate  the  variance-covariance  matrices  of  (1 
(P^)  and  of  2  ( Pj; )  from  S. 

Form  the  SEP  variance  approximation  using  equation 

(2)  . 

Computation  of  SEP  given  (i  and  E“  is  based  on  the  PL-1 
program  RAP  provided  by  JKU-APL.  RAP  is  a  general  program 
which  provides  the  radius  of  coverage  for  any  probability  in 
the  interval  (0,1)  for  either  2  or  3  dimensions.  Those 
portions  of  the  program  applicable  to  this  problem,  ie. 
probability  =  0.5  and  3  dimensions,  have  been  used  in  the 
FORTRAN  subroutine  FINDSEP  to  obtain  SEP.  Derivation  of  the 
equation  used  to  compute  SEP  is  given  by  L.  S.  Simpkins  [Ref. 
6]  and  is  implemented  in  the  subroutine  FINDSEP  using  Gaussian 
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Quadrature  (via  the  IMSL  subroutine  QTWODQ)  to  evaluate  the 
required  double  integral. 

Gaussian  Quadrature  is  also  used  to  evaluate  each  of  the 
double  integrals  in  the  W,  CC,  CS,  SC  and  SS  functions,  which 
are  required  in  the  evaluation  of  the  matrices  DTI,  DT2  and 
DT3. 


B.  TESTING 

The  purpose  of  deri  ing  the  second  order  approximation  of 
the  variance  of  SEP  is  to  yield  more  accurate  confidence 
intervals  for  SEP.  The  approximation  of  these  confidence 
intervals  is  derived  from  the  fact  that  since  (x  and  2  are 
maximum  likelihood  estimators  of  p  and  S  and  thus  SEP  is 
asymptotically  Normal  with  mean  SEP.  Thus,  for  example,  an 
approximate  95%  confidence  interval  for  SEP  is  given  by 


S^Ptl  .9eslvai{SEP) . 

In  order  to  test  the  affect  of  using  the  second  order 

approximation  for  the  variance  of  SEP,  tests  of  confidence 

interval  coverage  were  conducted  as  follows: 

-  Sample  3  0  random  vectors  from  and  compute 

p  and  2  for  the  sample. 

Calculate  SEP  from  ^  and  2  . 

Calculate  d^gp  by  replacing  jl  and  2  for  p.  and2 
in  both  the  first  order  and  second  order 
approximations  for  the  variance  of  SEP. 

Compute  the  large  sample  approximation  95% 

confidence  intervals  for  SEP  using  SEP±1.96^^. 
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Repeat  1000  times  to  calculate  the  percentage  of 
times  that  confidence  intervals  based  on  each 
approximation  cover  the  true  SEP. 

For  each  input  variance-covariance  matrix,  S, 
repeat  the  test  at  different  magnitudes  of  the 
mean. 

Trivariate  normal  random  vectors  are  obtained  using  the 
IMSL  subroutines  DCHFAC  and  DRNMVN.  DCHFAC  performs  Cholesky 
factorization  of  the  distribution  SIGMA  which  is  then  used  as 
an  input  to  DRNMVN  which  returns  the  desired  random  vector. 

The  estimates  of  and  Pj  are  found  by  replacing  ^  and 

S  with  (1  and  2  in  the  formulations  of  P^  and  P^  given  in 
Chapter  I.  The  results  of  this  test  are  summarized  in  Table 


TABLE  1.  95  PCT  CONFIDENCE  INTERVAL  COVERAGE 


SIGMA 

MU 

FIRST  ORDER 

SECOND  ORDER 

Cl 

Cl 

n 

0.8 

0.9 

0 

94.6 

95.2 

1.0 

0.8 

0 

mmm 

0.8 

1.0 

0 

1.0 

92.8 

93.3 

1.2 

0.8 

100 

-160 

270 

0 

93 . 3 

94.2 

-160 

400 

-480 

0 

270 

-480 

900 

0 

100 

95.6 

95.8 

250 

500 

10,000 

3 , 000 

2,000 

0 

93 . 3 

93.8 

3,000 

10,000 

3 , 000 

0 

2 , 000 

3 , 000 

10,000 

0 

8,000 

10,000 

12,000 

93 . 5 

94.0 
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These  results  indicate  that  there  is  little  accuracy  to  be 
gained  by  using  the  second  order  variance  approximation. 
There  is,  however,  a  situation  in  which  the  second  order 
approximation  becomes  important.  This  case  has  been  described 
by  K.  V.  Kitzman  [Ref.  3]  and  is  paraphrased  here. 

Weapons  system  targeting  can  be  measured  from  separate 
subsystems  (such  as  navigation,  guidance,  postboost,  etc.). 
Suppose  tnat  there  are  n  independent  subsystems  and  let 

,11)  represent  the  errors  for  the  subsystem 

i  =  1, . . . ,n,  then  the  impact  errors  for  the  weapons  system  are 

n  J,J.cI 

^  and  -  ^(^j,=n|ji,  2j.=nS) . 

1*1 

If  the  Y's  are  measured  directly  and  used  to  estimate Hj, 

2  „ 

and  2„  (to  get  the  system  SEP)  then  P„  =  — -  =  —  S  and 
^  m  m 

Py  =  — =  -^^(E®S).  If  on  the  other  hand  the  individual 

n 

X's  are  available  from  the  subsystems  then  where 

i  =  l 


J  =  1 


—  2  and 
m 


Then  P, 


O  n 

=  — (E0I1).  The  significant  factor  is  that  the  estimate  of 


differs  from  the  first  estimate  by  a  factor  of  n  while  the 

estimate  of  P^^  remains  unchanged.  Thus  as  the  number  of 

subsystems  becomes  larger  and  the  number  of  observations 
becomes  smaller  the  size  of  Pj;  can  become  very  small  in 

relation  to  P^, .  As  the  mean  approaches  0  it  follows  that 


dSEP 

a^l 


approaches  0  as  well  and  thus  the  contribution  of  P^  to 


the  estimate  can  become  insignificant  even  though  it  is  large 
in  relation  to  the  P^  term. 

This  affect  has  been  examined  as  follows: 

Divide  p  and  S  by  n,  the  number  of  iid  subsystems 
being  simulated. 

Draw  m  random  samples  from  the  distribution  given 
by  this  new  jji  and  S  for  each  of  the  n  subsystems. 
Calculate  p  and  £  for  each  of  the  n  subsystems. 
Sum  the  n  estimates  to  get  a  total  p  and  2 ,  the 
estimate  for  mean  and  variance  in  Y. 

Compute  confidence  intervals  for  SEP  using  both 
the  first  and  second  order  approximations  with 

P„  =  -E  and  P^  =  — (EOS). 


This  test  was  conducted  with  n  =  4  and  m  =  5.  The  results 
of  this  test  are  given  in  Table  2.  As  can  be  seen,  the  effect 
is  greatest  when  the  mean  is  0,  and  decreases  as  the  mean 
becomes  larger,  which  is  what  would  be  expected. 
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TABLE  2.  95  PCT  Cl  COVERAGE,  MULTIPLE  SUBSYSTEMS 


SIGMA 

MU 

FIRST  ORDER 

SECOND  ORDER 

Cl 

Cl 

100  -160  270 

0 

92.4 

97.6 

-160  400  -480 

0 

270  -480  900 

0 

5 

91.4 

94 . 8 

10 

20 

25 

89.4 

92.2 

50 

80 

50 

92 . 6 

94.4 

150 

200 

To  test  the  sensitivity  of  the  second  order  variance 

approximation  procedure  to  the  assumption  that  the  detonations 

are  normally  distributed,  tests  were  conducted  in  which  random 

vectors  were  drawn  from  a  mixed  normal  distribution.  The  test 

is  conducted  in  the  following  manner: 

Prior  to  each  sample  obtain  random  variable  p  from 
U(0,1)  distribution. 

If  p  is  less  than  .98  then  draw  a  sample  from 
and  proceed. 

if  p  is  greater  than  .98  then  draw  a  sample  from 
N{\i,  10*11)  and  proceed. 

Continue  until  30  samples  have  been  drawn,  then 
compute  all  values  and  proceed  as  above. 

This  procedure  simulates  a  distribution  with  the  same  mean 
but  larger  tails  than  the  normal  distribution.  The  confidence 
interval  coverages  are  then  compared  to  those  obtained  in  the 
base  case  from  the  same  input  distribution.  The  results  of 
this  test  are  summarized  in  Table  3. 
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TABLE  3 


95  PCT  Cl  COVERAGE,  MIXED  DISTRIBUTION 


IV.  CONCLUSIONS 


Examination  of  the  results  in  Table  1  shows  that  the 
second  order  approximation  procedure  provides  a  slight 
improvement  in  the  95%  confidence  interval  coverage.  The  real 
utility  of  this  method  is  brought  out  by  the  results  of  the 
multiple  subsystem  test  shown  in  Table  2,  where  confidence 
interval  coverage  improves  by  as  much  as  5%.  Since  this  case 
represents  the  actual  situation  in  missile  testing  (i.e.  very 
small  mean  and  independent  system  observations)  it  is 
reasonable  to  state  that  use  of  this  procedure  will  improve 
the  estimate  of  the  variance  in  SEP.  Furthermore  it  seems 
highly  unlikely  that  there  would  be  any  advantage  to  be  gained 
by  pursuing  a  higher  order  estimate  of  the  variance  of  SEP. 

The  mixed  distribution  test  shows  that  the  second  order 
approximation  is  only  slightly  more  robust  than  the  first 
order  approximation  if  the  distribution  is  actually  from  a 
mixed  and  not  a  true  normal.  It  also  indicates  that  if  the 
variance  is  extremely  large  then  the  affect  of  the  mixing  is 
greater  than  for  the  smaller  values  and  that  only  a  few 
extreme  values  are  required  to  greatly  reduce  the  confidence 
interval  for  the  estimate. 
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APPENDIX  A 


In  order  to  derive  the  equation  for  the  second  order 
approximation  for  variance  of  SEP  (equation  (1))  we  start  by 
notinq  that  the  second  order  Taylor  Series  Expression  for  the 
SEP  function  is  given  by  (ignoring  the  remainder  term): 

where  X,  ^SEP,  and  D^  j^SEP  are  defined  in  the  text,  X  is  the 

vector  consisting  of  the  maximum  likelihood  estimators  for 
each  element  of  X  and  AX  =  X-X .  Our  goal  then  is  to  find  an 

approximation  for  the  variance  of  SEP(X).  Because  SEP(X)  is 
asymptotically  unbiased  and  the  mean  square  error 

e[{SEP{X)-SEP{X))^]  =  Var(SEF{X))HE[SEP{X)]-SEP(X)f , 

we  can  get  this  approximation  by  looking  at  ELSITPCX) -SEP(X)  ]  ^ . 
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Applying  the  tensor  identity  ABC  =  A®C'^B.  to  the  last  term 

and  then  transposing  we  can  write  the  second  order  Taylor 
equation  as: 


SEP{X)-SEI\\)  ~  + 


"(AXOAX) 


Now  we  note  that  since  £(X)  =  1.  it  follows  that  £i;AA,)  =  0  . 
If  we  square  each  side  and  take  the  expected  value  we  see 
that : 

£[EEP(X)-EEi^A.)]^  =  [D^  ■^uSEP(X)]'^E[AXAX'^[D^  ^^SEi\k)] 

ruSEi^X)]  ^£[(A  A.® A  A.)(A  A,® A  A.)  „u5Ei^X)] 

In  order  to  evaluate  £({AX®AA.)(AX®AA)’]  we  note  that 

var(AX®AA)  =  2  (SR)  (P®P)  (SR)^  [Ref.  7:p.  44]  where  the  matrix 

SR  has  the  following  qualities  [Ref.  7:  pp.  32-38]: 

SR  is  symmetric 
-  (SR) (SR)  =  SR 

2  ( SR )  (tJ®y)  =  ij^Y  *  Y^U 
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for  example,  if  u  and  v  are  of  dimension  2  then 


SR 


Then  since, 

£[(AX0AA,)(AX<S>AX)T  =  var(AX0AX)+£(AX0AX)[£:(AX0AX)]^ 

and 

AXAX^  =  AX0AX 


if  we  define  P  to  be  £[AXAX^  =  var(AX)  we  see  that 

rar(AX0AX)+£i:AX0AX)[^AX®AX)]’'  =  2SR(i^P){SR)^->-££l 

=  2(SR)'^(F0P){SK)*E£l. 
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since  is  a  symmetric  matrix  it  follows  that 


(SR)D^,  ^SEP  =  dI^xSEP  and  we  see  that; 

^SEF{X)-SEP{X)f  »  [D^^j,SEP{X)YPID^^^SEP{X)] 


which  in  turn  is  approximately  equal  to  the  variance  of 
SEP( X )  . 

Finally,  since  2  is  a  3x3  symmetric  matrix,  we  need  only 
describe  the  uniqtae  terms  (which  we  have  denoted  by  S“)  to 
arrive  at  the  .--_cond  order  equation  for  SEP  given  by  equation 
(1)  . 
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APPENDIX  B 


A.  APPROACH 

It  is  seen  that  each  of  the  12  functions  consists  of 

a  double  integral  of  cos  and  sin  functions  multiplied  by  the 
function  g(i?sin0cos4>#  i?sin0sin(j>,  RcosS) .  Each  of  the  required 
derivatives  carries  through  the  integrals,  so  that  we  need  to 
evaluate  each  in  respect  to  the  function  g. 

We  will  first  evaluate  each  of  the  partials  of  g,  and  then 
carry  these  through  to  the  functions  . 

Recall , 

gi(i?,0,4>)  =  exp|--i[(i?sin0cos(J)-n^)2o^^ 

+  2(Rsin0cos<t»-|ii)(i?sin0sin(()-ji2)o^^ 

+  2(i?sin0cos4)-Hi)(i?cos0-H3)o^' 

+  (Rsin0sin<|>-H2)^®^^ 

+  2(Rsin0sin<t)-ji2)(^cos0-vi3)CJ^^ 

+  (Rcos0-}i3)‘o^^j} . 
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If  we  define  the  vectors  A^,  A^,  A3  and  A^  by 

sin^6cos^<j) 
2sin^6cos4)sin4) 
2sin0cos0cos(t) 
sin^0sin^<J) 
2sin0cos0sin<{) 

COS-0 

sin0cos<l)^i3 

sin0sin<{)ji3+sin0cos4)^2 
cos0|i3+sin0cos4)^3 
sin0sin4)|j.2 
cos0^i2  +  sin6si^<J)^l3 
cos0p.3 

A3  =  [iil  2^l3^l2  2|i3|i3  \il  2|i2^3  ^i3]^ 

i?sin0cos4>-|i3 
A^  =  i?sin0sin<})-[i2 
i?cos0-H3 
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Then  treating  R  as  an  independent  parameter  and  taking  the 
partial  with  respect  to  R  we  have: 

=  -sr-|a’-^(i?sin^0cos^4)-|jijsin6cos4)) 

+o^^(2Rsin^0cos<t»sin<t)-^jSin6sin4>-H2sin6cos4)) 

+o^^(2i?sin0cos0cos4>-HiCos0-ti.3sin0cos<t)) 

+  o^^(Rsin^0Gin^4>-ji2sin0sin(j)) 

+  o^^(2Rsin0sin<{)cos0-ji2Cos0-H3sin0sin(|)) 
+o^-(Rcos^0-ji3Cos0)| 


similarly,  we  derive: 


=  g-(o“/Rsin0cos4)-^3)+a^^(Rsin0sin4)-^2)'^o^^(-^cos0-^3)} 


dg 

a^2 


=  g-|a^^(Rsin0cos<t)-|j.3)+a^^(/!’sin0sin4>-pi2)'^o^^(Rcos0-p.3)} 


dg 

aiij 


g-|o^^(Rsin0cos4)-H3)  +  o^^(i\sin0sin(l)-ji2)‘^cr^^(^cos9-H3)} 
g-{[a” 
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In  order  to  derive  the  partials  with  respect  to  sigma  we 


will  first  generate  the  partials 


da 


and  the  partials  of 


kl 


with  respect  to  .  These  values  are  given  below, 


do 


11 


d 

da 


12 


da 


13 


da 


22 


do 


23 


da 


33 


sll 

o22 

-o33ol3 

^33  -o22oll 

g^ 

~®33_2o12o12 

-203^032 

-2o22a32 

g^ 

g' 

-2a“o“ 

-2o,3  ,, 

-2a2^o 

g^ 

g" 

g" 

<^33  _ollo22 

-o12<,22 

-o“o22 

-022a22 

g^ 

_2oiio 

23  _2<,12o23 

_2o”o23 

-2o22o23 

g^ 

g^ 

g^ 

*^22  __11_33 

“®12  _rtl2o33 

—  n 

-^-o22o33 

-  u  u 

g^ 

g^ 

w  V 

g^ 

iiii 

fill 

a 

”‘^23.  _„23oll 

^22_p33oll 

g' 

_oii 

ao^i 

■'  '■■  V/  VI 

g^ 

2g 

a 

^13_Oo23o12 

^^^32  _2o33o12  . 

-a22 

aoi2 

g^ 

g^ 

g 

a 

^32  -2o23o” 

g^ 

-2o33o33 

-033 

aoi3 

g 

a 

-o23o22 

11_q33q22 

-o22 

a022 

g^ 

2g 

a 

'''“-2a23o23 

g^ 

-2o33o23 

-023 

ao23 

g 

a 

__23-33 

-0^^0^3 

-a33 

aa33 

VI  VI 

2q 
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We  can  now  use  these  values  to  find  the  desired  partials 

=  g^-||^[(i?sin0cos<|)-pi)2(a^^)2 

+2(i^sin0cos<^)-^iJ^)(i^sin0sin<|>-^l2)o^^a^^ 

■t-2(i?sin0cos4>-jii)(;?cos0-p3)a’-^o^’- 

+(i?sin0sin<|)-ji2)^|o^^o^^“-^Y  j 

+2(i?s  in0sin4>-H2)(-^c:os0 -10,3)^02^0^^ +  -^^j 

+(i?cos0-H3)2^o23o“--^'  --^1 


=  g-|-|[o^ii?2(E  -1)  “A3  -2i?oii(2  -^)  “^2  ■')  “^^3] 

-  ?||[o“(E-r*(o  0  0^:^ 


Following  this  same  methodology  it  is  easy  to  see  that; 


p-  =  gjo^2(2-i)>vYo 

^Oi2  I  \  2g2 


P-  =  o22(2-i)“.i 


<^33 

*^23  Q  *^13 

2g2 

2g2  2g2  g 

-<^23 

®22  *^13  “*^12 

2g2 

2g2  g2  2g2 

-O33 

n  n  n  ^11 

KJ  U  \J  ' 

Q  g2 

[R^A,-2RA^^A,]-^ 
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023(S-^) 


^A^-2RA2+A^  - 


dg 

00,3 


\ 

0  0 

y 


[R^A^-2RA2+A^]- 


We  are  now  able  to  calculate  each  of  the  partials  of  thefv^^ 

by  simply  computing  the  effects  of  the  four  A  vectors  on  the 
original  W  functions  and  applying  the  formulas  derived  on  the 
previous  pages.  We  will  carry  through  the  complete  derivation 
for  and  then  show  the  solutions  for  the  other  cases.  Since 

the  vector  A3  does  not  involve  trigonometric  functions,  it 
will  not  be  necessary  to  compute  any  partials  for  it. 
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B.  PARTIALS  OF  W1 


= 


1  1  1  ~  f  f  cos0cos4>g(i?, 6, <J))d0d(j)  so  we  can  see  that 

'  '  '  J  Q  J  0 


the  integrals  of  gA  will  yield: 


cos0cos<l)g-Ai 


sin^0cos0cos^<J> 

2sin^0cos0cos^<()sin4) 

2sin0cos^0cos^<j) 

sin^0cos0cos<|)sin^<j) 

2sin0cos^0sin<|)cos<t) 

cos^0cos<|> 


^(^■^1, 1, 1, 1  1, 1,  1, 3, 1  ■'■^■^3, 1, 3,  l) 

,1,2,1  , 1 , 2 ,  l) 

^^1 , 1 , 1 , 1  “  ^^3 ,1,1,1'  ^^1 , 1 , 3 , 1  ^^3 ,1,3,1 
■^■^l ,  1 , 1 . 2  ~  ■^■^3 , 1 , 1 , 2 
^^3, 1,1,1 


f'’cosdcos<t>g-A^ 
J  0  j  0 


sin0cos0cos2(t>|ii 

sin0cos0cos<t>sin<l>^^+sin0cos0cos^(j>^2 

cos^0cos(J)iij^+sin0cos0cos^4)p.3 

sin0cos0cos<|)sin<J>^2 

cos^0cos<j)ji2+sin0cos0cos<J)sin<})^3 

COS^0COS<J)H3 


2Hi5'C3  2,2,1 
1  •^•^1 , 2 , 1 , 2 ^  2  ■^^i ,  2 , 2 , 1 
4^iC’C2  11  i'''2jl3SC3  2,2,1 
^^2‘^'^l,2,l,2 

4  |A2^^2 , 1 , 1 , 1  ^^3'^'^1 ,2,1,2 

^^3^^2, 1,1,1 
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i?sin0cos6cos^<j)-cos0cos<j)[ij^ 

i?sin0cos0cos(t)sin4)-cos0cos4)|i2 

i?COS^0COS4>-COS0COS<|)fl3 

^sin20cos^<j>-li3cos0cos4) 

p 

—  sin20sin24>-H2Cos0cos4> 
i?cos^0cos4>-n3cos0cos4> 

^  , 2 , 2 , 1 ^  1  , 1 , 1 , 1 

•^•^1,2,1,2'^~^^^2^^1, 1, 1, 1 
4  CC2 , 1, 1,  CC3 ,1,1,1 


1 

4 


f  "  Pcos0cos<|)g-A4  =  f  'fg 

J  0  j  0  J  Q  J  0 


For  future  use,  we  point  out  that  the  partials  of  A2  andA^ 

involve  only  3  distinct  trigonometric  terms,  which  we  will 
define  by; 


B  = 


Si 

■S2 

S3 


=  sin0cos4) 
=  sin0sin4) 
=  COS0 


then  we  see  that 


Sl^l 

S2^1+S3^2 

S3^1+S3^3 

^2^^2 

S3^2^S2|l3 

S3P3 
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and 


^4 


Thus,  we  need  only  compute  the  effects  of  and  B  on  each 
of  the  12  W  functions  to  obtain  the  desired  partials. 

C.  PARTIALS  OF  W2 


/'  2k  /•  K 

/  cos36cos4)Sf(R, 6, <|))d0d(})  so  the  partials  are 
0  Jo 


given  by; 


cos30cos4)Aj^ 


sin^6cos36cos^<t> 

2sin^0cos30cos^4>sin<|) 

2sin0cos0cos30cos^(l) 

sin^0cos30sin^(|)cos<|> 

2sin0cos0cos30sin<j>cos(f> 

COS"0COS30COS4> 


2CC. 


1,3, 3,1 


- rc  - rc 

'-'“1,5, 3,1  '“'-1,1,3,: 


2 ( 2  CS3 ^  2.  i.i~ *“ "^i ,5,1,1“ ^"^i ,  1 , 1 , 1  “ 2 ,  3 , 3 , 1  ,5,3,3  +  CS3 ,1,3,1) 

^(•^^1,5, 2,1  “^^1,1,2, 1) 

^  *“  ^1 ,3,1,1“  ^^1 ,5,1,1“  ^^1 ,1,1,1“^  ^^1 ,3,3,1'*'  ^^1 ,5,3,1'*’  ^^1 ,1,3,1 

,5,1,2  “^-^i  ,1,1,2 
^^1,5, 1,1  '*2  ^^1, 3, 1, 1  '*’^^1, 1, 1, 1 
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sin0cos30cos^4> 

cos30cos4>B  =  sin0cos30cos4>sin4> 
COS0COS30COS(J) 


1 

4 


^(^^1 ,4,2, 1  ,2,2 

■^■^1,4, 1,2  “•^■^1,2,1, 
2(CCi ,  4 , 1 , 1  CC^  2,1 


D.  PARTIALS  OF  W3 


^3  ~  CS-^  ,1,1,1  cos0sin4>gt-^^,  0,  <}))id0d<J)  so  the  partials  are 

given  by: 

sin^0cos0cos^4)sin<i) 
2sin^0cos0cos(j>sin^<l> 
2sin0cos^0cos(J>sin(J> 
sin^0cos0sin^<}> 

2sin0cos^0sin^<|) 
cos^0sin<|) 


cos0sin4)^j^  = 


C'^l.l, 

,l.l-^^3. 

1,1,1  “^"^i,  1, 

3 , 1  ^"^3 ,1,3,1 

2(cq., 

,  1 . 1  ~  ^^1 

,1,3,1  ”^^3,1 

,  1 , 1  ^^3 ,1,3,1 

SSi 

1,1, 2  “'^■^3,1, 

1,2 

CSi, 

l,3,l“^'^3,l. 

3,1 

2(SSi 

,1,2,1~^'^3,1 

,2,l) 

^■^3 .1,1,1 

cos0sin4)fl  = 

sin0cos0co3<|>sin(|) 

sin6cos0sin^4) 

1 

J 

•^•^1 ,2,1,2 
^  ■^•^1 ,2,2,1 

cos^03in<t) 

** 

^■^2, 1, 1, 1, 
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E. 


PARTIALS  OF  W4 


cs 


1,3, 1,1 


cos30sin<l)gf(i?,  6 ,  <J))d8d4) 


so  the  partials 


are  given  by: 


cos36sin4»Aj^ 


sin^0cos30cos^(|>sin<l> 

2sin^0cos30cos<t>sin^(t) 

2sinOcos0cos30cos<J>sin4) 

sin^0cos30sin^<l> 

2sin0cos0cos30sin^<J) 

cos~0cos30sin4) 


1,3,3, 

1  ~  ^"^l  ,5,1,1 

*  , 5 , 3 , 1  ■ 

-CS. 

X  / 

'1,3,3 

I'CCi  5  1 

1  ^^1 ,5,3,1 

-cq 

•^•^1,5,1,2~ 

‘^‘^1,1, 1,2 

2CSj 

.,3,3,l“^'^l. 

,3,1 

^(•^•^1, 5,2,1' 

■■^■^l,l,2,l) 

5, 1,  l'*'2C’'^l, 

3,1,1'*’  ^'^1 , 1 

,1,1 

cos30sin<J)B 


sin0cos30cos4)sin4)' 

sin0cos30sin^<|) 

cos0cos30sin<j> 


1 

4 


•^•^1 , 4 , 1 , 2  “•^•^1 , 2 , 1 , 2 
^(■^■^1,4,2,  l“‘^‘^l,2,2,  l) 
.^(^•^1 , 4  ,  1 ,  1  ^-^1 , 2 , 1 ,  l). 
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F.  PARTIALS  OF  W5 


=  SC,  10  0^  P*  r’'sin0si(i?,  6,  (|>)ded<j)  so  the  partials 

'  '  '  J  Q  J  0 


given  by: 


sin^0cos^4) 

2  s  in^0c  os<|)  s  i  n<l> 
2sin^0cos0cos<i> 
sin’0sin^<}> 
2sin^0cos0sin<J) 
cos-0sin0 


SC3 

1,2,1 

SS3 

1,1,2 

2(cq 

,1,1,1 

~^^3, 1,1,1 

SS3_ 

1,2,1 

2(CSi 

,1,1,1 

“  ^"^3 ,1,1,1 

SCi_ 

1,0.0 

"■^^3, 1,0,0 

sin^0cos4> 

2  ^^2 ,1,1,1 

sin0S  = 

sin^0sin4> 
sin0cos0 . 

1 

2 

2  ^^2 ,1,1,1 
^^1 ,2,0,0  . 

are 
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G 


PARTIALS  OF  W6 


Wg  =  SC,  2  0  0  ^  f  f  3in2dg(R,6,(^)ddd4>  so  the  partials  are 
given  by: 

sin^6sin26cos^<j) 

2sin^6sin26cos<J)sin<i) 
2sin6sin20cos0cos4> 

sin20Ai  = 

sin^esin2esin^4> 

2sin0sin20cos0sin4) 
cos^0sin20 

^  ,2,2,1“  ,4,2,1 

2  •^•^1 , 2 , 1 , 2  ~  •^■^1 , 4 , 1 , 2 
^(^^1 ,2,1, 1  ~^^1, 4 ,1,1  '*■2  SC2 ,1,1,1) 

2 ‘^•^1 , 2 , 2 , 1  “ ‘^•^1 , 4 , 2 , 1 
^(^•^1 ,2,1,1“  ^“^1 ,4,1,1  ■'■2  -552 ,1,1,1) 

^  "^^l  ,2,0,0  ■'■■S’C^^  ,4,0,0 

^^1 ,1,1,1“  ^^1 ,3,1,1 
—  1,1  '^‘^1,3,1,! 

■^^1 , 3 , 0 , 0  ■*■  "^^i  ,1,0,0 


sin0sin20cos4> 
sin20S  =  sin0sin2flsin4> 
sin20cos0 
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H.  PARTIALS  OF  W7 


=  SC^  3  0  0  =  f  f  sin30g(i?, 0, <}))ded(J)  so  the  partials 
'  '  '  j  0  j  0 


given  by: 


sin30Aj^ 


sin^0sin30cos^<|> 

2sin^0sin30cos(J>sin(t> 

2sin0sin30cos0cos<l> 

sin^0sin30sin^4) 

2sin0sin30cos0sin<J) 

cos^0sin30 


2  ,3,2,1“  ,1,2,1“  "^^l  ,5,2,1 

2  ■^•^1 , 3 , 1 , 2  “ ‘^•^1 , 1 , 1 , 2  “  •^■^1 , 5 , 1 , 2 
1  ^(^^l,l,l,l”^^l,5,l,l) 

4  2SSj^  3  2,i“5'Sj^  1  2,i“*^-^1,5,2,1 

^(^■^1 , 1 , 1 , 1  “ ^"^l ,  5, 1 ,  l) 
■^^1,5,0,o''’^'^^1,3,0,o'*’'^^1,1,0,0 


sin30S 


sin0sin30cos<t> 

sin0sin30sin4) 

sin30cos0 


2 


^^1 ,2,1,1 

^•^1,2,1,! 


-cc 

'-'-1,4, 1,1 
“  ^"^l  ,4,1,1 
'*’‘^^1, 2, 0, 0, 


are 
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I. 


PARTIALS  OF  W8 


Wo 


SC. 


2, 1, 1, 1 


sin^0cos4>g(i?,  6, 4>)d0d4> 

Jo  Jo 


so  the  partials  are 


given  by: 


sin^6cos<t>Jlj^ 


sin^0cos^<j> 

2sin^0cos^<i)sin<j> 

2sin^0cos0cos^<t) 

sin^0sin^<j>cos(i) 

2sin'0cos0cos(j)sin4) 

cos^0sin-0cos(t> 


1 

8 


1,3,1 

1 6  -  SS4  3 

2(250^  2,2,1~’^^1,4,2,i) 
^(•^^4 , 1 , 1 , 1  “■^^4 , 1 , 3 . 1) 
2  •^•^1 ,2,1,2  ,4,1,2 

8(SC2  1  1  1  1) 


sin^0cos4)S 


sin^0cos^<|) 

sin^0sin<^cos<j) 

sin^0cos0cos4> 


^  ,1,2,1 

■^■^3 ,1,1,2 

(^^1 , 1 , 1 , 1  ~  ^^3 ,1,1,1) 
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J.  PARTIALS  OF  W9 


■^^3  1  2  1  ^  f  f  sin^0cos^4)g(J?,  6,<|))d0d4>  so  the  partials  are 
'  '  '  J  0  J  0 


given  by: 


sin^0cos^4)A3^  = 


sin^0cos‘<j) 

2sin-0cos^<J>sin<|> 

2sin^0cos0cos^(}> 

sin-0sin^<l)cos^<{> 

2sin‘0cos0cos^<|)sin({) 

cos-0sin^0cos“<i) 


1,4,1 

^■^5 , 1 , 1 , 4  ■*"  2 ,  1 , 1 , 2 

®  ( ^^1 , 1 , 3 , 1  ^*"3 ,1,3,1''’  ^^5 ,1,3,1) 

^(■^•^5 , 1 , 2 , 1  "•^•^5 , 1 , 4 , 1) 

^(^■^1, 1,  1,  1  “^"^l,  1, 3,1  “2  C'Sj  J+2CS3  3^  J  +  CS5  3^  3_ 

^(•^^3 , 1 , 2 , 1  “"^^5 , 1,  2 , 1) 


sin‘‘0cos^<l)  1,3,1 

sin^0cos''<t)5  =  sin^0sin4)cos^<|>  =  ^(•^■^4.i,i.i“'^-^4, 1,3,1) 

sin’0cos0cos^<J), 
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K. 


PARTIALS  OF  WIO 


/  r  ^ 

/  sin^0sin4>g(i?,6,4>)(:i0d<t>  so  *cne  partials  are 
0  Jo 

given  by: 

sin^6cos^<t>sin4) 

2sin^0cos(j)sin^<i) 

.  .  2sin^0cos0cos(f)sin(b 

sin20sin4)/^.  = 

sin^0sin^4) 

2sin^0cos0sin“4) 
cos“0sin“0sin<J) 


^(■^■^4 , 1 , 1 , 1  “ , 1 , 3 ,  l) 

1 6(50^ ,  1  _  1 , 1  ~SC^ ,  1 , 3 ,  i) 
8  8SS,,,,3,, 

^  •^•^1 , 2 , 2 , 1  ”  2  •^•^1 ,4,2,1 

8(S52_  3  1  -SS^  3  1,  l) 


sin^0cos4)sin4) 

■^■^3 ,1,1,2 

sin^0sin4)5  = 

sin^0sin^<J)  =  — 

2  ‘^•^3  ,1,2,1 

sin^8cos0sin<})J 

^(^■^1,1,  1,  ,  1, 1,  l). 
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L. 


PARTI ALS  OF  Wll 


=  SS^  112“/'  f  sin^0sin24)gi(i?,  6, 4»)d0d4)  so  the  partials 
are  given  by: 

sin^6cos^<J)sin2<J) 
2sin^0cos(i)sin<j)sin24t 
.  .  ,  2sin‘0cos0cos4)sin2(J) 

sin^6sin24>Aj^  = 

sin®0sin  '})sin2(l) 
2sin^0cos0sin<l)sin24) 
cos^0sin^0sin2<J) 


■^•^5, 1, 1,4  1, 2 

2(2SC5j^2,i“-^^5,1,1,4“‘^^5,1,1,2) 

,  1 , 1 , 1  “  2  ^"^3 , 1 , 1 , 3  “  2  ^“^3 , 1 , 1 , 1 
2  ,  1 , 1 , 2  “  •^*^5 , 1 , 1 , 4 

,  1 , 1 , 3  “  2  ^^3 ,1,1,i'*'2^^3,1,1,3''’^^5 
‘^(*^■^3 , 1 , 1 , 2  “'^‘^5 , 1 , 1 , 2) 

sin^0cos4>sin24>  ‘^(■^■^4, 1,1, 3 '^‘^■^4, 1,1,1) 

sin^0sin2<j)S  =  sin^0sin4>sin2<J>  = 

sin^0cos0sin24)  '^^^x,2,x,2~^^x,i.x.2 


,1,1,3"**  ^*^5 ,1^1,1) 
,  1,1,1  ~  ^^5, 1,1, 3) 


1  ^(‘^■^l,l,l,3‘'‘‘--^l 

4 

^(^^i,i,i,i“^^i 
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M.  PARTIALS  OF  W12 


W..  = 


'12 


/’  2k  /*  X 

/  sin^0sin^<j>p(i?,  0,<|))d0d<J>  so  the  partials 
0  Jo 


are  given  by: 


sin^0sin2<t)i^j^ 


sin®0cos^<|)sin^4) 

2sin®0cos(j)sin^<i> 

2sin^0cos0cos4)sin^(J) 

sin®0sin‘(|) 

2sin^0cos0sin^<|) 

cos^0sin^0sin^(|) 


■^(■^•^5 , 1 , 2 . 1  “■^'^5 , 1 , 4  ,  l) 

2 ‘^■^5 , 1 , 1 , 2  ~  ■^‘^S ,  1 , 1 , 4 

®  , 1 , 1 , 1  “  , 1 , 3 , 1  ~  2  CC3 , 1 , 1 , 1  ■*■  2  CC3 ,1,3^3  +  CCg  _  1  ^  1  ^  1  “  CC5 , 1 , 3 , 1 ) 

3  4  3 

,1,3,1  “2  C53  ^3,3^3  +  CSg  ,3,3^3) 

^  ( ■^‘^3 , 1 , 2 , 1  ” '^■^5 , 1 , 2 , 1 ) 


sin^0cos4>sin^4> 

®('^Q,i,i,i“'^Q,  1,3,1) 

sin^0sin^<t)S  = 

sin^0sin^<}> 

1 

8 

^  '^■^4 ,1,3,1 

sin^0cos0sin^<l). 

.  2 ‘^‘^1,2, 2,1  "■^■^1,4, 2,1  , 
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APPENDIX  C 


^  =  ^6 
5*^5  2wl 

du^  =  _J_ 

8Wg  2  ^5 


and  all  other  partials 


0. 


PARTIALS  OF  U4,..,U9  WRT  MU 


2Vc 


2  Vc 


11°  o22 


*^6  22 
-x - 

2Wc 


0^3 

SUg 

d\i 

50 


2 


Partials  of  U5 


du 


da 


—  = 


11 


du, 


5  _ 


da 


=  -2a^^u. 


12„  _  1 


12 


5  ? 


’33 


\  ”s 


'23 


[  4^5 


IV'.; 


du, 


5  _ 


3o 


■2a'^^u^-- 


13 


i  4lV^5  Vs 


V.;  V.; 


8u 


^  =  -a^^u  - 


da 


22 


5  2 


i?Vl,  _  tx,v,o 

2  V.-  v. 


8u, 


5  _ 


da 


23 


-2o“u5-^ 


'23 


^Vji -2^3  v^Q  \ 


V.; 


-a 


13 


^»^12~t^2t*^lo'l_g  i  _  1^3^10 


Vc 


12! 


4  V.. 


V. 


5  } 


du^ 

da 


-a^^u^-- 


33 


'12| 


^^12  f^2^10 


V.:  V. 


-o 


22 


5  / 


2  V.;  v. 


3.  Partials  of  U6 


du, 


da 


*  =  -o“u. 


11 


-2a^"u, 


4v5  2V5 


4  Vs  2V5  ^ 


52 


53 


^^7  22 


^  .  -2o»> 


diLj 

~ 


^  <7^  \  ^**^5  25^5  j  4^5  2^5  j 

'  gn  4  2w,  j  2^3  2k^3  j 


5.  Partials  of  U8 


^  =  -o^^u.- 


4  Wj  2  ^5 


4  W5  2  W5 


=  -20^2^1 


^  fR(w^-w2)_ 


®  g2|  4^5  2^5  J  '“I  4^5  2^5 


.2o33^^-^oJ:^^r\)-iii^ 
®  q^  “I  2w^3 


4  V5  2  V5 


4(^5  2(v^3 


—  = 


=  -2o^^Ug--^  o,,  -0.3 

0023  g'  I  4^3  2^5  I 


4w^5  2^5 
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?  7  — 


— r  a 


4W.,  2W.3 


4w^  2Wc  4W.;  2w^ 


6.  Partials  of  U9 


—  =  O23 

/T* 


8W5  4^5 


8  V5  4  V5 


.20^2^  -J-O 

®  12[  4;^^  2ur3 


R(W^-W^)_  [  ^'^l-^2)_  ^1^6 


8  V5  4  V5 


8  Vg  4  V5 


^^L  -  -2013.  _ J_  o  -o 

'i3  "  4W5  J  8V5  4^5 


i2(V3-V4)  lljVg 


=  -o^“Uq- — -a 


^  g2  13|^  0^^  j  11|^  0^^  4,^,^ 


u.-X[o..('^-'>-‘-'>>-!^Vo,i^‘^‘-'*'^>-i^ 


“g  2  11 


8  (^5  4  V5 


8V5  41^5 


^  =  -o”u. 
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D.  PARTIALS  OF  U4 


»  •  •  » 


U9  WRT  R 


3tl 

E.  PARTIALS  OF  U4,..,U9  WRT  W 
1.  Partials  of  U4 

du^  ^  i?a33 
dw^  8^5 

dW2  8  ^5 


56 


1 


2.  Partials  of  U5 


8U5 


8u4 

8h^5 


U.. 


^^4  ^  Ra^^ 

dWg  2  V5 


du^  = 

4^5 


^^5  ^  Ra^^ 

dw^  4  V5 

^^5  ^  _  j?0^^ 

8^4  4 


=  -A„, 

8W5  W^  ■ 


^^5  ^  J?0^^ 

8wjj^  21^5 
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J?o^2 


du^ 


3.  Pairtials  of  U6 

^ 

dWj^  4  V5 

^  =  ^£211 

dw2  4  V5 

=  i?o^^ 

8^3  4(^5 

=  -£sll 

dw^  4  V5 


8^5 


4^5 


4 ' 

8^  = 

8W7  4  Vj 
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4 


Partials  of  U7 


0^3  8  W5 


^ 

dw^  8  V5 


dw^  V5 


0^33  4^5 

8^12  2(^5 


5.  Partials  of  U8 

8^3  4  V5 


^ 

8w'2  4  h^5 


^^8  _  Ro^^ 
81*^3  4  V5 
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dUg  ^  _Ra^ 

dw^  4  V5 


dw^  V5  ®  4  V5 


4  ■  -^(<'"•‘.-“^•^*””>‘3) 


dug  = 
dw^  4  V5 


6.  Partials  of  U9 

du^  = 

3wj,  8V5 


^^9  ^ 

dw^  8  V5 


8^3  8  h^5 


du^  = 

8(v'^  8  w^ 


^  =  -_iu  ^  ^q” 
8^5  *S  ^  ®  'S 
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dUg 


8Wj  8  ^5 


^"3) 
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APPENDIX  D 


PROGRAM  SEPCMP 

*  PROGRAMMER:  LT  ARTHUR  F.  BROCK 

*  DATE:  12  MAY  1991 

*  LAST  MODIFIED:  29  AUGUST  1991 

*  PROGRAM  DETERMINES  FIRST  AND  SECOND  ORDER  APPROXIMATIONS  TO 

*  VARIANCE  OF  SEP  GIVEN  VALUES  OF  MU  AND  SIGMA  (VARIANCE-COVARIANCE 

*  MATRIX)  AS  INPUTS.  OTHER  INPUTS  ARE  PARAMETERS  FOR  USE  IN  THE 

*  IMSL  SUBROUTINE  QTWODQ  USED  TO  PERFORM  GAUSSIAN  QUADRATURE 

*  EVALUATION  OF  DOUBLE  INTEGRALS. 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

INCLUDE  'PMUCOM  DEF' 

INTEGER  I ,OT1 ,OT2,SAMPSZ,K 

REAL  DT1(10,9),DT2(22,10),DT3(9,22),SEPM(9.9),PMU(3,3>,PSIG(6,6) 

C  ,SEPMU(3),SEPSIG(6),CI(2),MPT<3,3) 

CALL  INIT(SAMPS2,MPT) 

URITE(16,51)  R.NUMTR 

WRITE(16,52)(MU(I), 1=1,3) 

DO  20  1=1,3 

WRITE(16,53)(MSIG(I,J),J=1,3) 

20  CONTINUE 
OT1=0 
OT2=0 

DO  10  I=1,SAMPS2 
CALL  INIT2(MPT) 

CALL  UVALS 

CALL  DT1VAL(0T1,SEPMU,SEPSIG) 

CALL  0T2VAL(DT2) 

CALL  DT3VAL(DT3) 

CALL  SEPEST(DT1,DT2,DT3,SEPM) 

CALL  PMUS(MSIG,PMU,PSIG) 

CALL  BN0EST(SEPM,PMU,PSIG,SEPMU,SEPS1G,CI) 

IF(SEPSET  .LT.  (R  CKD)  .OR.  SEPSET  .GT.  (R+CI(1)))  THEN 
OT1=OT1+1 

IFCSEPSET  .LT.{R-CI(2)).OR.  SEPSET  .GT . (R+CI (2) ) )  THEN 
OT2=OT2+1 
ENDIF 
ENDIF 

PRINT*, 'I  =  ',I,REAL(I-OT1)/REAL<I),REAL{I-OT2)/REAL(I) 
WRITE(16,50)  RCI(1),R+CI(1),OT1,R-CI(2),R+CI(2),OT2 
10  CONTINUE 

WRITE(16,55)  SAMPSZ,REAL(OT1)/REAL(SAMPSZ).REAL{OT2)/REAL(SAMPSZ) 

50  FORMAT(1X,2( '  ( ' , F7.2, ' , ' , F7.2, ' ) ' ,2X, I3,2X)) 

51  FORMATdX, 'TEST  FOR  SEP  =  ', F7. 2,/, IX, ' SAMPLES  PER  TRIAL  =  ',12) 

52  FORMATdX,'  MU  ',/ ,3(  F7. 2, ','),/,/,  IX, 

C  '  SIGMA') 

53  FORMATdX, 3(F9. 3, 3X)) 

55  FORMATdX,/, 'SAMPLE  SIZE:  ' ,  IA,3X,2(F5.3,8X)) 

STOP 

END 
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SUBROUTINE  INIT<SAMPSZ,MPT ) 


INITIALIZE  FILES,  READ  INPUT  DATA  AND  INITIALIZE  VARIABLES. 

INCLUDE  'COM  DEF' 

INCLUDE  'PMUCOM  DEF' 

INTEGER  N.SAMPSZ.RK 

REAL  MPT(3,3),TOL,NUMSIG(3,3) 

EXTERNAL  RNSET.DCHFAC 
OPEN(15,FILE=  /MUSIG  DATA  AT) 

0PEN(16,FILE='/0UTDT  DATA  AT) 

READ(15,»)  MU(1) 

READ(15,*)  MU(2) 

READ(15,*)  MU(3) 

READdS,*)  SIGd) 

READdS,*)  SIG<2) 

READdS,*)  SIG{3) 

READdS,*)  SIG(4) 

READdS,*)  SIGCS) 

READdS,*)  SIG{6) 

READdS,*)  ERRABS 
READdS,*)  ERRREL 
READdS,*)  IRULE 
READdS,*)  NUHTR 
READdS,*)  SAMPSZ 
READdS,*)  ISEED 
READdS,*)  NUMSIM 

CALL  RNSET(ISEED  +  INT(HUd  )  +  SIGd  ))) 

DET=SIGd)*(SIG(4)*SIG(6)-(SIG(S)**2))-SIG(2)*(SIG(2)*SIG(6) 

C  ■SIG(3)*SIG<S))+SIG(3)*(SIG<2)*SIG(5)SIG(3)*SIG{4)) 

IF  (OET  ..E.  0.0000000001)  THEN 
URITEd6,S0)  DET 
STOP 
END  IF 

SIGINVd)=<SIG(4)*SIG(6)-(SIG(S)**2))/0ET 

SIGINV(2)s(SIG<3)*SIG<S)-SIG(2)*SIG(6))/OET 

SIGINV(3)»(SIG(2)*S1G<S)-SJG(3)*SIG(4))/DET 

SIGINV<4)  =  (SIGd)*SIG(6)-(SIG(3)**2))/DET 

SIGINV(S)s<SIG(3)*SIG(2)-SIGd)*SIG(S))/0ET 

SIGINV(6)  =  (SIGd)*SIG(4)-{SIG(2)**2))/DET 

DENOM=<SORT{2.0*PI)**3)*SQRT<DET) 

MSIGd,1)=SIGd) 

MSIGd,2)=SIG(2) 

MSIGd,3)=SIG(3) 

HSIG(2,1)=SIG(2) 

MSIG(2,2)=SIG(4) 

MSIG(2,3)=SIG(S) 

MSIG(3,1)=SIG(3) 

MSIG(3,2)=SIG{5) 

MSIG(3,3)=SIG(6) 

CALL  FINDSEP 

MSIGd,1)=SIGd)/NUMSIM 

MSIGd  ,2)=SIG(2)/NUMSIM 

MS!Gd,3)=SIG(3)/NUMSIM 

MSIG(2,1)=SIG(2)/NUMSIM 

MSIG(2,2)=SIG(4)/NUMSIM 

MSIG(2,3)=SIG{S)/NUMSIM 

MSIG<3,1)=SIG(3)/NUMSIM 

MSIG(3,2)=SIG{S)/NLIMSIM 

MSIG(3,3)=SIC;6)/NUMSIM 

TOL=100.0*DMACH(4) 

CALL  0CHFAC(3,MSIG,3,T0L,RX,MPT,3) 

CALL  SETDAT 
CALL  MAKES 

50  FORMATdX, 'ILLEGAL  INPUT  MATRIX,  DETERMINANT  =  ',F8.5,'  PROGRAM', 
C  '  TERMINATED.') 

RETURN 

END 
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SUBROUTINE  INIT2(MPT) 


INITIALIZE  VALUES  USED  EACH  TIME  THROUGH  LOOP. 

REAL  MPT(3,3),TMU(3),TSIG(6) 

INCLUDE  'COM  DEF- 
INCLUDE  'PMUCOM  DEF* 

20  CONTINUE 

SIG(1)=0.0 

SIG(2)=0.0 

SIG(3)=0.0 

SIG(4)=0.0 

SIG(5)=0.0 

SIG(6)=0.0 

MU(1)=0.0 

MU(2)=0.0 

MU(3)=0.0 

DO  15  I=1,INT(NU«SIM) 

CALL  TRIALS(NUNTR,MPT,TMU,TSIG) 

SIG(1)=SIG(1)+TSIG(1) 

SIG(2)=SIG(2)+TSIG<2) 

SIG(3)=SIG(3)+TSIG<3) 

SIG(4)=S;G(4)+TSIG(4) 

SIG(5)=SIG(5)+TSIG(5) 

SIG(6)=SIG<6)+TSIG(6) 

MU(1)=MU(1)+TMU(1) 

MU(2)=MU<2)+TMU(2) 

MU(3)=HU(3)+TMU(3) 

15  CONTINUE 

MSIG(1,1)=SIG(1) 

HSIG(1,2)=SIG(2) 

MSIG(1,3)=SIG(3) 

MSIG(2,1)=SIG(2) 

MSIG(2,2)=SIG(4) 

MSIG(2,3)=SIG(5) 

MSIG(3,1)=SIG(3) 

MSIG(3,2)=SIG(5) 

MSIG(3,3)=SIG(6) 

DETsSIG(1)*(SIG(4)*SIG(6)-<SIG(5)**2))-SIG(2)*(SIG(2)*SIG(6) 
C  -SIG(3)*SIG(5))+SIG(3)*(SIG<2)*SIG(5)-SIG{3)»SIG(4)) 

IF  (DET  .GT.  0.00000000001)  GOTO  30 
GOTO  20 
30  CONTINUE 

SIGINV(1)={SIG(4)*SIG(6)<SIG(5)**2))/DET 

SIGINV(2)=<SIG<3)*SIG(5)-SIG(2)*SIG<6))/OET 

SIGINV(3)={SIG<2)*SIG(5)SIG(3)*SIG<4))/OET 

SIGINV<4)=<SIG<1)*SIG<6){SIG(3)**2))/DET 

SIGINV(5)*(SIG(3)*SIG<2)SIG<1)*SIG(5))/OET 

SIGINV{6)={SIG(1)*SIG(4)-(SIG<2)»*2))/0ET 

DENOM=(SQRT{2.0*PI)**3)*SQRT(OET) 

CALL  FINDSEP 

RETURN 

END 


SUBROUTINE  MXTMLT(M1 ,M2, ANS, I , J,K) 

PERFORMS  MATRIX  MULTIPLICATION  Ml < I , J)»M2( J ,K)  TO  PRODUCE 
OUTPUT  MATRIX  ANS(I,X) 

INTEGER  L,ROW,COL,I,J,K 
REAL  M1(I,J),M2(J,K),ANS(I,X) 

DO  10  ROW=1,I 
DO  20  COL=1,X 

ANS(ROW,COL)=0.0 
20  CONTINUE 
10  CONTINUE 

DO  30  ROW=1,I 
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DO  40  C0L=1,K 
DO  50  L=1,J 

ANS(ROU,COL)=ANS(ROW,COL)+H1(ROU,L)*M2(L,COL) 
SO  CONTINUE 

40  CONTINUE 
30  CONTINUE 
RETURN 
END 


SUBROUTINE  SEPESTTDTI ,DT2,0T3,SEPM) 

PRODUCES  THE  9  BY  9  MATRIX  OF  PARTIALS  OF  SEP 
INCLUDE  'COM  DEF' 

REAL  DT1{10,9),DT2(22,10),DT3(9,22),SEPM(9,9).TEMP(22.9) 
CALL  MXTMLT(DT2,DT1,TEMP,22,10,9) 

CALL  MXTMLT(DT3,TEMP,SEPH,9,22,9) 

RETURN 

END 


SUBROUTINE  SETDAT 

SETS  THE  INPUT  VALUES  OF  MU,  SIG  AND  SEP  TO  BE  HELD  CONSTANT 
THROUGHOUT  THE  LOOPS  OF  THE  PROGRAM. 

INCLUDE  'COM  DEF' 

INCLUDE  'PMUCOM  DEF> 

INTEGER  I,J 
D010  1=1,3 

MUSET(I)=MU(I) 

10  CONTINUE 
SEPSET=R 
RETURN 
END 


SUBROUTINE  TRIALS(N,MPT,TMU,TSIG) 

DRAWS  N  RANDOM  SAMPLES  FROM  NORMAL (MU, SIGMA)  AND  DETERMINES 
ESTIMATES  FOR  MU  AND  SIGMA  BASED  ON  THESE  SAMPLES. 

INCLUDE  'COM  DEF* 

INCLUDE  'PMUCOM  DEF' 

INTEGER  I,N,J 

REAL  MPT(3,3),RVAR(50,3),RSQ(50,3),RMLT(50,3),RHN(6),OP(3) 

C  ,NUMFAC,TSIG(6),TMU(3) 

EXTERNAL  ORNMVN , RNSET , DMACH , DRNUNF 
DO  10  1=1, N 

CALL  DRNMVN(1,3,MPT,3,OP,1) 

RVAR( I , 1 )=OP( 1 )+MUSET( 1 )/NUMSIM 

RVAR(I,2)=OP(2)+MUSET(2)/NUMSIH 

RVAR(I,3)=OP(3)+MUSET(3)/NUMSIM 

RSO{I,1)=RVAR(I,1)**2 

RSO(I,2)=RVAR(I,2)**2 

RSO(I,3)=RVAR(I,3)**2 

RMLT( I , 1 )=RVAR{ I , 1 )*RVAR( 1,2) 

RMLT( I ,2)=RVAR{ I , 1 )»RVAR( I ,3) 

RMLT( I ,3)=RVAR( I ,2)*RVAR( I ,3) 

10  CONTINUE 

TMU(1)=FMN(RVAR,1,N) 

TMU(2)=FMN(RVAR,2,N) 

TMU(3)=FMN{RVAR,3,N) 

RMN<1)=FMN(RSQ,1,N) 

RMN(2)=FMN(RS0,2,N) 

RMN(3)=FMN(RSQ,3,N) 

RMN(4)=FMN(RMLT,1,N) 
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RMN(5)=FMN(RMLT,2,N) 

RHN(6)=FMN(RHLT,3,N) 

NUHFAC=1.0 

TS1G{1)=(RHN(1)-TMU(1)**2)*NUMFAC 

TSIG(2)=(RHM(A)-TMU(1)*TMU(2))*MUMFAC 

TSIG{3)=(RMN(5)-TMU(1)*TMU(3))*MUHFAC 

TSIG{4)=(RHN(2)-TMU(2)**2)*NUMFAC 

TSIG(5)=(RMN{6)-TMU(2)*TMU(3))*NUMFAC 

TSIG(6)=(RMN(3)-TMU{3)**2)*NUMFAC 

RETURN 

END 


REAL  FUNCTION  FHN(VEC,COL,N) 

DETERMINES  THE  MEAN  OF  N  ITEMS  IN  COLUMN  COL  OF  VECTOR  VEC 

INTEGER  I.N.COL 
REAL  VEC<50,3),T0T 
TOT=0.0 
DO  10  1=1, N 

TOT=TOT+VEC(I,COL) 

10  CONTINUE 

FMN=TOT/REAL(N) 

RETURN 

END 


REAL  FUNCTION  G(X) 

REAL  X 

G=0.0 

RETURN 

END 


REAL  FUNCTION  H(X) 

REAL  X 

H=6. 2831854 

RETURN 

END 


SUBROUTINE  WVALS 

DETERMINES  VALUES  OF  U,  CC,  CS,  SC  AND  SS  FUNCTIONS  USING  THE 
IMSL  SUBROUTINE  OTWODO. 

REAL  G,H,H1,H2,ERREST 
INCLUDE  'COM  OEF' 

INCLUDE  'WVEC  DEF' 

EXTERNAL  G, H,DTWOOO,CC1 1 1 1 ,CC131 1 , CC31 1 1 , CS1 1 1 1 ,CS131 1 , CS31 1 1 , 

C  SC1100,SC1200,SC1300,SC2111.SC3121,SS2111,SS3112,SS3121, 

C  CC1113,CC1131,CC1211,CC1331,CC1411,CC1511,CC1531,CC2111, 

C  CC3113,CC3131,CC5111,CC5113,CC5131,CS1113,CS1131,CS1211, 

C  CS1331,CS1411,CS1511,CS1531,CS2111.CS3113,CS3131,CS5111, 

C  CS5113,CS5131,SC1121,SC1212,SC1221,SC1321,SC1400,SC1412, 

C  SC1421,SC1500,SC1521,SC3100,SC4111,SC4113,SC4':31,SC5112, 

C  SC5114,SC5121,SC5141,SS1112,SS1121.SS1212,SS1221,SS1312, 

C  SS1321,SSl412,SS1421,SS1512,SS1521,Sb4111,SS4113,SS4131, 

C  SS5112,SS5114,SS5121,SS5141 

B=PI 

CALL  DTWOOQ  (CC1 1 1 1 , A,B , G, H.ERRABS , ERRREL , IRULE ,W( 1 ) , ERREST ) 

CALL  DTUOOO  (CC1311 , A,B,G,H,ERRABS,ERRREL, IRULE,W(2),ERREST) 
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CALL  DTWC»Q  (CS1 1 1 1 , A,B, G, H,ERRABS,ERRREL, IRULE ,W(3), ERREST ) 
CALL  DTWOOQ  {CS131 1 , A,B,G, H.ERRABS, ERRREL, IRULE,W(4),ERREST ) 
CALL  DTUOOQ  (SC1100, A, B,G, H.ERRABS, ERRREL, IRULE, W(S), ERREST) 
CALL  DTWOOQ  (SC1200, A, B,G, H.ERRABS, ERRREL. IRULE, W(6), ERREST) 
CALL  DTWOOQ  (SC1300, A, B, G, H.ERRABS, ERRREL, IRULE, W{7) .ERREST ) 
CALL  DTWOOQ  (SC21 1 1 , A, B,G, H.ERRABS, ERRREL, IRULE, W(8), ERREST ) 
CALL  DTWOOQ  (SC3121 , A, B,G, H.ERRABS, ERRREL, IRULE, W(9), ERREST) 
CALL  DTWOOQ  (SS21 11 , A, B,G, H, ERRABS, ERRREL, IRULE, W(10), ERREST ) 
CALL  DTWOOQ  (SS31 12. A, B, G, H.ERRABS, ERRREL. IRULE. W{1 1 ), ERREST ) 
CALL  DTWOOQ  (SS3121 , A, B.G, H.ERRABS, ERRREL, IRULE, W(12), ERREST) 
CALL  DTWOOQ  (CC1 1 13, A, B.G, H.ERRABS, ERRREL, IRULE. CC(1 ), ERREST) 
CALL  DTWOOQ  (CC1 131 , A, B.G, H.ERRABS, ERRREL, IRULE, CC(2), ERREST) 
CALL  DTWOOQ  (CC121 1 , A, B, G, H.ERRABS. ERRREL, IRULE ,CC(3) .ERREST ) 
CALL  DTWOOQ  (CC1331 , A, B, G, H.ERRABS, ERRREL, IRULE. CC(5), ERREST ) 
CALL  DTWOOQ  (CCU11 , A, B.G, H.ERRABS, ERRREL. IRULE. CC{6), ERREST) 
CALL  DTWOOQ  (CC1511 , A, B.G, H.ERRABS, ERRREL. IRULE, CC{8). ERREST) 
CALL  DTWOOQ  (CC1531 , A, B.G, H.ERRABS, ERRREL. IRULE, CC(9), ERREST) 
CALL  DTUOOQ  (CC21 1 1 ,A,B,G,H, ERRABS, ERRREL. IRULE, CC{10). ERREST) 
CALL  DTWOOQ  (CC31 1 1 , A, B.G, H.ERRABS, ERRREL, IRULE. CC{17), ERREST) 
CALL  DTWOOQ  (CC31 13, A,B, G, H.ERRABS, ERRREL, IRULE. CC(1 1 ) .ERREST ) 
CALL  DTWOOQ  (CC3131 , A, B.G, H.ERRABS, ERRREL, IRULE ,CC( 12), ERREST ) 
CALL  DTWOOQ  (CC51 1 1 , A, B.G, H , ERRABS, ERRREL, IRULE, CC( 14), ERREST ) 
CALL  DTWOOQ  (CC51 13, A, B , G, H.ERRABS, ERRREL , IRULE, CC( 15), ERREST ) 
CALL  OTWOOQ  (CC5131 , A, B, G, H.ERRABS, ERRREL, IRULE ,CC(16), ERREST ) 
CALL  DTUOOQ  (CS11 13, A, B.G, H.ERRABS, ERRREL. IRULE, CS(1 ), ERREST ) 
CALL  DTUOOQ  (CS1 131 , A, B.G, H.ERRABS, ERRREL. IRULE, CS(2), ERREST) 
CALL  DTUOOQ  (CS121 1 , A, B.G, H, ERRABS, ERRREL, IRULE, CS(3) .ERREST ) 
CALL  DTUOOQ  (CS1331 , A, B.G, H.ERRABS, ERRREL, IRULE, CS(5). ERREST) 
CALL  DTUOOQ  (CS141 1 , A, B, G, H .ERRABS, ERRREL , IRULE ,CS(6), ERREST ) 
CALL  OTWOOQ  <CS151 1 , A, B.G, H.ERRABS, ERRREL , IRULE ,CS(8), ERREST ) 
CALL  OTWOOQ  (CS1531 , A, B ,G, H.ERRABS, ERRREL, IRULE. CS(9). ERREST ) 
CALL  OTWOOQ  <CS21 1 1 , A, B.G, H.ERRABS, ERRREL, IRULE, CS( 10), ERREST ) 
CALL  DTWOOQ  (CS31 1 1 , A, B.G, H, ERRABS, ERRREL. IRULE, CS( 17), ERREST ) 
CALL  DTWOOQ  (CS31 13, A, B.G, H.ERRABS, ERRREL, IRULE, CS(18), ERREST) 
CALL  DTUOOQ  (CS3131 , A, B.G, H.ERRABS, ERRREL, IRULE. CS(11 ), ERREST) 
CALL  DTUOOQ  (CS51 1 1 , A, B.G, H.ERRABS, ERRREL, IRULE. CS( 13), ERREST) 
CALL  DTUOOQ  (CS51 13, A, B.G, H.ERRABS, ERRREL, IRULE, CS(14) .ERREST ) 
CALL  DTWOOQ  (CS5131 , A, B.G, H.ERRABS, ERRREL. IRULE, CS(16), ERREST) 
CALL  DTUOOQ  (SC1 121 , A, B.G, H.ERRABS, ERRREL, IRULE, SC( 1 ), ERREST) 
CALL  DTUOOQ  <SC1212, A, B.G, H.ERRABS, ERRREL, IRULE, SC(18), ERREST ) 
CALL  DTUOOQ  (SC1221 , A, B.G, H.ERRABS, ERRREL. IRULE, SC(2), ERREST) 
CALL  OTWOOQ  (SC1321 , A, B.G, H.ERRABS, ERRREL, IRULE, SC(3), ERREST) 
CALL  OTWOOQ  (SC1400, A, B, G, H.ERRABS, ERRREL , IRULE,SC(4),ERREST ) 
CALL  OTWOOQ  (SC1412, A, B.G, H.ERRABS, ERRREL , IRULE, SC( 19) , ERREST ) 
CALL  DTWOOQ  (SC1421 , A, B.G, H.ERRABS. ERRREL , IRULE, SC(5), ERREST ) 
CALL  OTWOOQ  (SC1500, A, B, G, H.ERRABS, ERRREL . IRULE ,SC(6), ERREST ) 
CALL  OTWOOQ  (SC1521 , A, B.G, H.ERRABS, ERRREL , IRULE, SC(7), ERREST ) 
CALL  OTWOOQ  (SC3100, A, B.G, H.ERRABS, ERRREL , IRULE, SC(9), ERREST ) 
CALL  OTWOOQ  (SC41 1 1 , A, B.G, H.ERRABS, ERRREL , IRULE ,SC( 10), ERREST ) 
CALL  DTUOOQ  (SC41 13, A,B, G, H,ERRABS,ERRREL , IRULE,SC( 1 1 ) .ERREST ) 
CALL  OTWOOQ  (SC4131 , A, B, G, H .ERRABS .ERRREL , IRULE,SC( 12) , ERREST ) 
CALL  DTUOOQ  (SCS1 12, A, B, G, H.ERRABS, ERRREL, IRULE, SC( 13) .ERREST ) 
CALL  DTUOOQ  (SCS1 14 , A, B, G, H, ERRABS, ERRREL, IRULE, SC( 14) .ERREST ) 
CALL  OTWOOQ  (SC5121 , A, B, G, H .ERRABS, ERRREL, IRULE ,SC( 15), ERREST ) 
CALL  DTWOOQ  (SC5141 , A, B, G, H, ERRABS, ERRREL , IRULE ,SC( 16) , ERREST ) 
CALL  DTWOOQ  (SS1 1 12, A,B,G, H, ERRABS.ERRREL , IRULE ,SS( 1 ) .ERREST ) 
CALL  OTWOOQ  (SS1 121 , A, B.G, H.ERRABS, ERRREL , IRULE,SS(2) .ERREST ) 
CALL  DTUOOQ  (SS 1212, A, B.G, H.ERRABS, ERRREL , IRULE. SS(3) .ERREST ) 
CALL  OTWOOQ  (SS1221 , A, B.G, H.ERRABS, ERRREL, IRULE, SS(4), ERREST) 
CALL  DTUOOQ  (SS1312, A ,B ,G, H.ERRABS, ERRREL , IRULE.SS(5 ) .ERREST ) 
CALL  DTUOOQ  (SS1321 , A, B.G, H.ERRABS .ERRREL , IRULE,SS(6) .ERREST ) 
CALL  DTUOOQ  (SS1412, A, B.G, H , ERRABS.ERRREL , IRULE, SS(7), ERREST ) 
CALL  DTWOOQ  (SS1421 , A, B, G, H, ERRABS, ERRREL , IRULE , SS(8) .ERREST ) 
CALL  DTUOOQ  (SS1512, A ,B,G, H,ERRABS,ERRREL , IRULE,SS(9) , ERREST ) 
CALL  DTUOOQ  (SS1521 , A ,B, G, H.ERRABS .ERRREL , IROLE,SS( 10),ERREST ) 
CALL  DTWOOQ  (SS41 1 1 ,A,B , G, H, ERRABS .ERRREL , IRULE , SS( 12) .ERREST ) 
CALL  DTWOOQ  (SS41 13, A, B, G, H .ERRABS , ERRREL , IRULE , SS( 13) .ERREST ) 
CALL  DTUOOQ  (SS4131 , A, B, G, H .ERRABS , ERRREL. IRULE , SS( 14 ), ERREST ) 
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CALL  DTUOOQ  (SS5112,A,B,G,H,ERRABS,ERRREL,IRULE.SS(15),ERREST) 
CALL  DTUOOQ  (SS51 U, A, B,G, H,ERRABS,ERRREL, IRULE,SS(16),ERREST ) 
CALL  DTUOOQ  (SS5121 ,A,B.G, H,ERRABS,ERRREL, IRULE.SS(17),ERREST) 
CALL  DTUOOQ  (SS5U1 , A, B, G,  H,ERRABS,ERRREL,  IRULE,SS(  18) , ERREST ) 
RETURN 
END 


SUBROUTINE  FINDSEP 

BASED  ON  PL-1  PROGRAM  RAP,  PRODUCES  SEP  BASED  ON  HU  AND  SIGMA 
FINDSEP  DIAGONALIZES  SIGMA  BASED  ON  EIGENVALUES  AND  THEN  CALLS 
SUBROUTINE  ITERATION,  UHICH  DETERMINES  SEP 

INCLUDE  'COM  DEF' 

REAL  EIGVLS(3),EIGMAT(3,3),Q1{3,3),NMEAN{3),NEUSIG(3,3) 

C  ,TEMP(3,3),EIGNML<3,3) 

INTEGER  I,J 
EXTERNAL  DEVCSF 

CALL  DEVCSF(3,MSIG,3,EIGVLS,EIGMAT,3) 

CALL  ORfHOCEIGMAT.EIGNML) 

DO  11  1=1,3 
DO  21  J=1,3 

01(I,J)=EIGNML(J,I) 

21  CONTINUE 
11  CONTINUE 

CALL  MXTMLT(Q1,MU,NMEAN,3,3,1) 

CALL  MXTMLT(MSIG,EIGNML,TEMP,3,3,3) 

CALL  MXTMLT(01,TEMP,NEUSIG,3,3,3) 

CALL  ITERATION(NMEAN,NEUSIG) 

RETURN 

END 


SUBROUTINE  ORTHO  <A1,A2) 

ORTHONORHALI2ES  MATRIX  OF  EIGENVECTORS  OF  SIGMA 

REAL  A1(3,3),A2(3,3),L,TEMP,TEMP2 
INTEGER  I 

L=SQRT((A1(1,1)**2)*<A1<2,1)**2)+<A1(3,1)**2)) 

DO  10  1=1,3 

A2(I,1)=A1<I,1)/L 
10  CONTINUE 

TEMP=A2( 1 , 1 )*A1 (1 ,2)+A2(2, 1 )*A1 (2,2)+A2(3, 1 )»A1{3,2) 
DO  20  1=1,3 

A2(I,2)=A1(I,2)TEMP*A2(I,1) 

20  CONTINUE 

L=SQRT(A2(1,2)**2«A2(2,2)»*2+A2(3,2)**2) 

DO  30  1=1,3 

A2(I,2)=A2(I,2)/L 
30  CONTINUE 

TEMP=A2<1,1)»A1(1,3)+A2(2,1)*A1(2,3)+A2(3,1)»A1(3,3) 
TEMP2=A2(1,2)*A1<1,3)+A2(2,2)»A1(2,3)+A2(3,2)»A1(3,3) 
DO  40  1=1,3 

A2(I,3)=A1(I,3)-TEMP*A2(I,1)TEMP2*A2<I,2) 

40  CONTINUE 

L=S0RT(A2(1,3)*»2+A2(2,3)**2+A2(3,3)**2) 

DO  50  1=1,3 

A2(I,3)=A2(I,3)/L 
50  CONTINUE 
RETURN 
END 


SUBROUTINE  ITERATION(HEAN,SIGMA) 


*  GENERATES  SEP  BASED  ON  HU  AND  SIGHA  FOR  DIAGONALIZED  SIGHA 

INCLUDE  'COH  DEF' 

INCLUDE  'SEPCOM  DEF' 

REAL  TOL,MEAN(3),SIGMA(3,3),MOH1,TRACE(3),SQSIG(3,3).CUSIG(3,3) 
C  ,TEMP(3),C2,D0F,RL0W,TPL0W,RHIGH,TPHIGH.CRATI0, integral 

C  ,MOH2,MOH3 

TOL=0.001 

CXX»1.0/SIGHA(1,1) 

CYY=1.0/S1GMA(2,2) 

CZZ=1.0/SIGMA(3,3) 

CSANT=OUPITW«SQRT{CXX*CYY*CZZ) 

C=(CXX»(HEAN(1)**2)+CYY*(HEAN(2)»*2)+CZZ*(MEAN(3)»»2))/2.0 
CALL  MXTMLTCMEAN, MEAN, HOHl, 1,3,1) 

CALL  HXTMLT(SIGMA,SIGMA,SQS1G,3,3,3) 

CALL  MXTMLT(SIGMA,MEAN,TEMP,3,3,1) 

CALL  HXTMLT(MEAN, TEMP, M0M2, 1,3,1) 

CALL  MXTMLT(SIGMA,SQSIG,CUSIG,3,3,3) 

CALL  MXTMLT(SOSIG,MEAN,TEMP,3,3,1) 

CALL  MXTMLT(MEAN, TEMP, M0H3, 1,3,1) 

DO  10  1=1,3 
TRACE(I)=0.0 
MN(I)=MEAN(I) 

10  CONTINUE 
DO  20  1=1,3 

TRACE! 1 )=TRACE( 1 )+SI GMA( I , I ) 

TRACE(2)=TRACE(2)+SQSIG( 1,1) 

TRACE(3)=TRACE(3)+CUSIG( 1,1) 

20  CONTINUE 

M0M1=M0MUTRACE(1) 

M0M2=2 . 0*<  TRACE  <  2 )+2 . 0*MOM2 ) 

M0M3=8 . 0* ( TR ACE  <  3 )*3 . 0*MOM3 ) 

BETA=(M0M3**2)/(M0M2*»3) 

OOF=8.0/BETA 

C2=OOF»(1.0-<2.0/<9.0»0OF)))**3 

RAOIUS=SOfiT(A8S(SQRT{ABS(MOM2/(2.0*OOF)))*<C2-OOF)*MOM1)) 

RLOU=0.95*RAOIUS 

TPLO«=EVALT(RLOU) 

40  CONTINUE 

IF(TPLOW  .GT,  0.5)  THEN 
RLOW=0.9*RLOW 
TPLOW=EVALT(RLOW) 

GOTO  40 
END  IF 

RHIGH=1 .105»RLOW 
TPMIGH=EVALT(RMIGM) 

50  CONTINUE 

IF  (TPHIGH  .LT.  0.5)  THEN 
RHIGH=1 .1‘RHIGH 
TPHIGH=EVALT(RHIGH) 

GOTO  50 
END  IF 

CRATIO=<0.5TPLOW)/(TPHIGHTPLOU) 

RA0IUS=RL0W+(RHIGHRL0W)»CRATI0 

ICNT=0 

60  CCWTINUE 

ICNT=ICNT+1 

IF((RHIGH-RLOW)/RHIGH  .GT.  TOL)  THEN 
INTEGRAL=EVALT(RADIUS) 

IF(A8S( INTEGRAL  -  0.5)  .LT.  TOL)  GOTO  75 
IFdNTEGRAL  .GT.  0.5)  THEN 
RHIGH  =  RADIUS 
TPHIGH=INTEGRAL 
ELSE 

RLOW=RADIUS 
TPLOW= integral 
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END  IF 

RAO  I US=RLOW+(RH I GH - RLOW)*CRAT 10 
IF  (ICNT  .GE.  150)  THEN 
PRINT*, 'COUNT  EXCEEDED' 

GOTO  75 
END  IF 
GOTO  60 
END  IF 

75  CONTINUE 
R=RAOIUS 
RETURN 
END 


REAL  FUNCTION  EVALT(RAD) 

*  EVALUATES  THE  AREA  COVERAGE  BY  TRIVARIATE  NORMAL  FOR  INPUT  R 

INCLUDE  'COM  DEF' 

INCLUDE  'SEPCOM  DEF' 

REAL  ANS, LOWER, UPPER, RAD 

EXTERNAL  RAPFC3,RAPL0W,RAPUP 

RADIUS=RAO 

LOWER=0.0 

UPPER=6. 2831853 

CALL  OTWOOQ(RAPFC3,LOWER,UPPER,RAPLOU,RAPUP,ERRABS,ERRREL,IRULE, 
C  ANS,ERREST) 


EVALT=ANS 

RETURN 

END 


REAL  FUNCTION  RAPLOW(X) 

REAL  X 

RAPLOW=0.0 

RETURN 

END 


REAL  FUNCTION  RAPUP(X) 
REAL  X 

RAPUP=3. 1415927 

RETURN 

END 


REAL  FUNCTION  RAPFC3(THETA,PHI ) 

*  FUNCTION  USED  TO  FIND  SEP 

REAL  CP, SP, CT, ST, AA,SA, SAMI, B8,B2DA,ERFARG,EXPARG, THETA, PHI 
INCLUDE  'COM  DEF' 

INCLUDE  'SEPCOM  DEF' 

EXTERNAL  ERF 
CP=COS(PHI) 

CT=COS(THETA) 

SP=SIN(PHI) 

ST=SIN(THETA) 

AA=(CXX*SP*SP*CT*CT*CYY*SP*SP*ST*ST*CZZ*CP*CP)/2.0 

SA=SQRT(AA) 

SAH1=1.0/SA 

BB=(CXX*MN(1)*SP*CT+CYY*MN(2)*SP*ST+CZZ*MN(3)*CP)/2.0 

B2DA=BB*BB/AA 

ERFARG=SA*RA0IUS+SAM1*BB 

EXPARG= - AA*RAO I US*R AO  I  US  -  2 . 0*BB*R AO  I US • B2DA 

IF  (EXPARG  .GT.  -32.0)  THEN 

F=ERF{ERFARG)*SAM1*(0.5*B2OA)-EXP(EXPARG)*(RAOIUS-BB/AA) 
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C 

ELSE 

F=ERF{ERFARG)*SAM1*(0.5+B20A) 

ENDIF 

ERFARG=SAH1*BB 

EXPARG=-B20A 

IF  (EXPARG  .GT.  -32.0)  THEN 

F=F-ERF(ERFARG)*SAM1»(0.5+B2DA)-EXP(EXPARG)*(BB/AA)*SQPIINV 

ELSE 

F=F-ERF(ERFARG)*SAM1*(0.5+B20A) 

ENDIF 

EXPARG=B20A-C 

IF  (EXPARG  ,GT.  -32.0)  THEN 
F=F«CSANT*EXP(EXPARG)*SP/AA 
IF(F  .LT.  0.000000000000001)  F=0.0 
ELSE 
F=0.0 
ENDIF 
RAPFC3=F 
RETURN 
END 


SUBROUTINE  DT1VAL  (DT1 ,SEPMU,SEPSIG) 

SETS  THE  VALUES  OF  THE  MATRIX  DTI  AND  THE  VECTORS  SEPHU  AND 
SEPSIG. 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

REAL  DT1(10,9),SEPMU<3),SEPSIG(6) 

INTEGER  I,J 
DO  100  1=1,9 
DO  200  J»1,9 

IF  (I  .NE.  J)  THEN 
OT1<I,J)=0.0 
ELSE 

DT1(I,J)=1.0 

ENDIF 

200  CONTINUE 
100  CONTINUE 

DT1(10,1)=U(8)/U(5) 

0T1(10,2)=W{10)/U(5) 

0T1(10,3)*W(6)/<2.0*W<5)) 

OT1(10,A)*{R/(8.0*W(5)))*(A.O»SIGINV(1)*W(9)*2.0»SIGINV(2)*W(11) 

C  ♦SIGINV(3)*<W(1)-U<2)))-(W<8)/(2.0»W(5)))* 

C  (SIGINV(1)*MU(1)*SIGINV(2)*MU{2)+SIGINV(3)»MUC3)) 

DT1(10,5)«(R/(4.0*W<5)))*(2.0»SIGINV(1)*W(11)*4.0*SIGINV(2)*W(12) 
C  ♦SIGINV<3)*(W<3)-W(4)))-(U(10)/W(5))» 

C  (SIGINV<1)*MU{1)+SIGINV(2)»MU<2)+SIGINV(3)*MU(3)) 

DT1(10,6)={R/<4.0*W(5)))*(SIGINV(1)*(W(1)-W(2))*SIGINV(2)*(WC3) 

C  -W(4))-^SIGINV(3)*(W(7)*U(5)))-(W(6)/(2.0*W{5)))* 

C  (SIG1NV(1)»MU(1)*SIGINV(2)»MU<2)+SIGINV(3)*MU(3)) 

DT1(10,7)=(R/(8.0*W(5)))»(2.0*SIGINV(2)»W<11)*4.0*SIGINV(4) 

C  •W(12)»SIGINV(5)*{U(3)-W(4)))-<W(10)/(2.0*W(5)))* 

C  (SIGINV<2)»MU{1)+SIGIMV(4)»MU(2)*SIGINV(5)»MU(3)) 

DT1(10,8)=(R/(4.0»U(5)))*(SIGINV(2)*(W(1)-W(2))*SIGINV(4)*(W(3) 

C  -W(4))*SIGINV(5)*(W(7)*W{5)))-(U<6)/(2.0»W(5)))* 

C  (SIGINV(2)*MU(1)+SIGINV<4)*MU(2)+SIGINV(5)*MU(3)) 

OT1(1O,9)={R/(8.0*U(5)))*(SIGINV(3)*(W(1)-W(2))*SIGINV(5)*(W(3) 

C  -W(4))*SIGINV(6)*{U(7)*W(5)))-(W(6)/(4.0*W(5)))* 

C  (SIGIMV(3)*MU(1)+SIGINV(5)*MU(2)+SIGINV(6)*MU(3)) 

DO  10  1=1,3 

SEPMU(I)=OT1(10,I) 

10  CONTINUE 
DO  20  1=1,6 

SEPSIG(I)=DT1(10,I+3) 
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20  CONTINUE 
RETURN 
END 


SUBROUTINE  DT2VAL(0T2) 

SETS  THE  VALUES  OF  THE  MATRIX  DT2.  THE  VECTORS  A1  AND  B1  ARE 
IN  THIS  SUBROUTINE  AND  THEN  PASSED  TO  THE  SUBROUTINE  FIXHAT, 
WHICH  ACTUALLY  SETS  THE  VALUES  OF  0T2. 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

REAL  DT2(22,10),SIGM<3,3),TEMP(3),ANS(3),BU3),A1(6) 

INTEGER  I,J 
DO  100  1=1,10 
DO  200  J=1,10 

IF  (I  .NE.  J)  THEN 
DT2(I,J)=0.0 
ELSE 

0T2(I,J)=1.0 
END  IF 

200  CONTINUE 
100  CONTINUE 

SIGM{1,1)=SIGINV(1) 

SIGH(1,2)=SIGINV(2) 

SIGH(1,3)=SIGINV(3) 

SIGM(2.1)=SIGINV<2) 

SIGM(2,2)=SIGINV(4) 

SIGM(2,3)=SIGINV(5) 

SIGM(3,1)=SIGINV(3) 

SIGM(3,2)=SIGINV(S) 

SIGM(3,3)=SIGINV(6) 

B1(1)=SC<2)/2.0 

B1(2)=SS(3)/4.0 

B1<3)=CC(10) 

A1(1)=CC(2)CC(12) 

A1(2)=2.0»<W(3)-CS(17)-CS(2)+CS(11)) 

A1(3)=2.0*(SC<1)-W<9)) 

A1(4)=U(1)-CC(17)-CC(2)+CC(12) 

A1(5)=SS<1)-U(11) 

A1(6)=CC(17) 

CALL  FIXMAT(A1,B1,11,0T2,SIGM) 

B1(1)=(SC(5)-SC<2))/2.0 

B1{2)=(SS(7)SS<3))/4.0 

B1{3)=(CC{6)+CC(3))/2.0 

A1(1)=(2.0*CC(5)CC(9)-CC<2))/4.0 

A1(2)=(2.0*W{4)-CS(8)-W(3)-2.0*CS(5)+CS{9)+CS{2))/2.0 

A1(3)=(SC(7)SC{1))/2.0 

A1(4)=(2.0«W(2)-CC<8)-U<1)-2.0*CC(5)+CC(9)+CC(2))/4.0 
A1(5)=(SS(9)-SS(1))/4.0 
A1(6)=(CC(8)*2.0*W(2)*W<1))/4.0 
CALL  FIXMAT(A1,B1,12,DT2,SIGM) 

B1<1)=SS{3)/4.0 

B1(2)=SS(4)/2.0 

81(3)=CS(10) 

AU1)=W(3)-CS{17)-CS(2)*CS(11) 

A1{2)=2.0*(W{1)CC(2)-CC(17)+CC(12)) 

A1{3)=SS(1)U(11) 

A1(4)=CS(2)-CS(11) 

A1(5)=2.0*{SS(2)-W(12)) 

A1(6)=CS(17) 

CALL  FIXMAT(A1,B1,13,0T2,S1GM) 

B1(1)=(SS(7)SS{3))/4.0 

81(2)=(SS(8)SS(4))/2.0 

B1{3)=(CS(6)+CS(3))/2.0 

A1{1)=<2.P*W(4)-2.0*CS(5)CS(8)*CS(9)-W(3)*CS(2))/4.0 

A1(2)=(2.0*U(2)-2.0«CC{5)-CC(8)*CC(9)W(1)+CC(2>)/2.0 
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A1(3)=<SS(9)-SS(1))/4.0 
A1(4)=(2.0*CS(5)CS(9)-CS(2))/4.0 
A1(5)={SS(10)-SS{2))/2.0 
A1(6)=<CS(8)+2.0*U(4)+U(3))/4.0 
CALL  FIXMAT(A1,B1,14,DT2,S1GM) 

B1(1)=W(8) 

B1(2)=W<10) 

B1{3)=W(6)/2.0 

A1(1)=U(9) 

A1(2)=W(11) 

A1<3)=2.0*(U(1)-CC(17)) 

AU4)=U(12) 

A1(5)=2.0»(W(3>-CS(17>) 

A1(6)=U(5)SC(9) 

CALL  FIXMAT{A1,B1,15,DT2,SIGH) 

B1(1)»(W(1)-U(2))/2.0 
B1{2)*(U(3)-W(4))/2.0 
B1(3)*(U(7)+U{5))/2.0 
Al(1)-(2.0*SC(2)-SC(S))/4.0 
A1{2)»(2.0*SS<3)-SS(7))/4.0 
A1(3)=(2.0*U{8)+CC(3)-CC<6))/2.0 
A1(4)«(2.0*SS{4)-SS(8))/4.0 
Al(5)=(2.0*W(10)*CS(3)-CS<6))/2.0 
A1(6)=(2.0*W(6)+SC(4))/4.0 
CALL  FIXMAT(A1,B1,16,DT2,S1G«) 

B1<1)*(CC{3)CC(6))/2.0 
B1<2)=(CS(3)-CS(6))/2.0 
B1(3)*(SC<4)+U(6))/2.0 
A1(1)«(2.0*SC(3)SC<1)SC(7))/4.0 
A1(2)=<2.0*SS(5)-SS(1)SS(9))/4.0 
A1(3)*<U<1)CC(8))/2.0 
A1(4)=<2.0«SS(6)-SS(2)SS{10))/4.0 
A1(5)=(W(3)-CS(8))/2.0 
A1(6)=(SC(6)+2.0»W(7)*U(5))/4.0 
CALL  F1XMAT(A1,B1,17,DT2,SIGM) 

B1<1)*W(9) 

B1(2)»W{11)/2.0 

B1(3)»W{1)CC<17; 

A1<1)=SC(12) 

A1(2)=2.0*(SS(12)SS(14)) 

A1<3)=(2.0»SC(2)-SC<5))/4.0 

A1(4)=SC(10)SC(12) 

A1(5)=(2.0*SS(3)SS(7))/8.0 

A1(6)=U(8)-SC(10) 

CALL  FIXMAT(A1,B1,18,0T2,SIGM) 

B1(1)=SC(12) 

B1{2)=SS<12)-SS(14) 

B1{3)=(2.0*SC(2)SC(5))/8.0 

A1(1)=SC{16) 

A1(2)=(SS(16)*2.0*SS(15))/4.0 

A1{3)=2.0«{CC(2)-2.0»CC{12)+CC(16)) 

A1(4)=SS(17)SS(18) 

A1(5)=2.0»{U(3)CS(2)-2.0*CS(17)‘2.0*CS<11)+CS(13)CS(16)) 

A1(6)=W<9)-SC<15) 

CALL  FIX«AT(A1,B1,19,DT2,SIGM) 

B1{1)=W(11)/2.0 

B1{2)=W(12) 

B1(3)=W(3)-CS(17) 

A1{1)=SS(12)SS(U) 

A1(2)=2.0«(SC(10)-SC(12)) 

A1(3)=(2.0*SS(3)-SS(7))/8.0 

A1(4)=SS(14) 

A1<5)=(2.0*SS(4)-SS(8))/4.0 

A1(6)=W(10)SS(12) 

CALL  FIXMAT(A1,B1,20,DT2,SiGM) 

B1(1)=(SS(13)*SS(12))/2.0 
B1(2)=(SC(10)  SC(11))/2.0 
B1(3)=(2,0*SS(3)SS(7))/8.0 
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A1(1)=(SS{16)+2.0*SS(15))/4.0 
A1(2)=(2.0»SC{15)  'i:(H)  SC(13))/2.0 
A1(3)=CS(1)+U{3)-2.0*CS(18)-2.0»CS{17)+CS(U)+CS(13) 
A’(4)=(2.0*SS(15)-SS(16))/4.0 

A1{5)=W(1)-CC(1)-2.0»CC(17)+2.0*CC(11)+CC(14)-CC(15) 

A1(6)=W(11)-SS(15) 

CALL  FIXMAT(A1,B1,21,DT2,SIGH) 

B1{1)=SC<10)SC{12) 

B1(2)=SS(14) 

B1(3)=(2.0*SS(4)-SS(8))/a.0 

A1(1)=SS(17)-SS(18) 

A1(2)=(2.0»SS(15)-SS<16))/4.0 

A1(3)=2.0»(U(1)-CC(2)-2.0*CC(17)+2.0»CC(12)+CC(14)-CC(16)) 

A1(4)=SS{18) 

A1(5)=2.0*(CS{2)-2.0*CS(11)+CS(16)) 

A1(6)=W(12)SS(17) 

CALL  FIXMAT(A1,B1,22,0T2,SIGH) 

RETURN 

END 


SUBROUTINE  F1XMAT(A1 ,B1 ,N,DT2,SIGH) 

SETS  THE  VALUES  OF  0T2,  GIVEN  A1.  B1,  SIGHA  INVERSE  (SIGH) 
AND  THE  ROW  NIMBER  TO  BE  SET  AS  INPUTS. 

INCLUDE  'WVEC  DEF> 

INCLUDE  'COM  DEF' 

REAL  DT2(22,10),SIGM(3,3),A1(6),B1{3),A2(6),A3(6),A4(3),MSUM(6) 
C  ,G(6),ANS(3) 

INTEGER  N,I 
DO  100  1=1,3 

A4(I)=B1(I)*R-MU(I)*U<N-10) 

100  CONTINUE 

CALL  MXTMLT(SIGM,A4,ANS,3,3,1) 

DT2(N,1)=ANS(1) 

0T2(N,2)=ANS{2) 

DT2<N,3)=ANS{3) 

A2(1)=B1(1)*MU(1) 

A2(2)=B1<2)*MU(1)+B1(1)*MU(2) 

A2(3)=B1<3)*HU<1)+B1(1)»MU<3) 

A2(4)=B1(2)*HU(2) 

A2(5)=B1(3)*MU(2)+B1(2)*MU(3) 

A2(6)=B1{3)*HU<3) 

DO  200  1=1,6 

HSUM(I)=R*A1(I)A2(I) 

200  CONTINUE 

CALL  HXTHLT(SIGINV,HSUH,ANS, 1,6,1) 

DT2(N,10)=-ANS(1) 

A3(1)=W(N-10)»<MU<1)»»2) 

A3(2)=2.0*W(N- 10)*HU(1 )*MU(2) 

A3<3)=2.0*W(N-10)»MU(1)*HU(3) 

A3(4)=U(N-10)*(MU{2)**2) 

A3(5)=2.0*U(N-10)*MU(2)*HU(3) 

A3(6)=W(N- 10)*{MU(3)»*2) 

DO  300  1=1,6 

MSLIH<  I  )=(R**2)*A1  ( I )  •  2 .0*R*A2(  I  )+A3(  I ) 

300  CONTINUE 

G(1)=SIGINV(1)«SIGINV(1) 

G(2)=SIGINV{2)*SIGINV(1) 

G(3)=SIGINV(3)*SIGIHV(1) 

G<4)=SIGIHV{4)»SIGINV<1)-SIG(6)/DET 

G(5)=SIGIHV(5)*SIGINV<1)*SIG<5)/DET 

G<6)=SIGINV(6)*SIGINV(1)SIG(4)/0ET 

CALL  HXTMLT(G,HSUM,ANS, 1,6,1) 

0T2{N,4)=(ANS(1)-W(N-10)*SIGINV(1))/2.0 

G(1)=SIGIHV(1)»SIGIHV(2) 

G(2)=SIGIHV(2)*SIGIHV(2)*SIG(6)/(2.0*DET) 
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G(3)=S1G1NV(3)»SIGINV(2)-S1G(5)/(2.0*DET) 

G(4)=SIGINV(A)*SIGINV(2) 

G(5)=31G1HV(5)*SIGINV(2)-SIG(3)/(2.0*DET> 

G(6)=:SIG:HV(6)*SIGINV(2)*SIG(2)/DET 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,5)=ANS(1)-W(N-10)»SIGINV(2) 

G(1>=SIG1NV(1)*SIGINV(3) 

G(2)=SIGINV(2)*SIGINV(3)-S!G(5)/<2.0*OET) 

G(3)=SIGINV(3)*SIGINV(3)*SIG(4)/(2.0*0ET) 

G(4)=SIGIMV(4)»SIGINV(3)+SIG(3)/DET 

G(5)=SIGIMV(5)*SIGINV(3)-SIG(2)/(2.0*DET) 

G(6)=SIGINV(6)*SIGINV(3) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

DT2(N,6)=ANS(1)-W<N-10)»SIGINV(3) 

G(1)=SIGINV(1)*SIGINV(4)-SIG<6)/DEf 

G(2)=SIGINV(2)*SIGINV(4) 

G(3)=SIGINV(3)*SIGINV(4)+SIG{3)/DET 

G(4)=SIGINV(4)*SIGINV(4) 

G(5)=SIGINV(5)*SIGINV<4) 

G(6)=SIGIMV(6)*SIGINV<4)-SIG{1)/0ET 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

0T2(N,7)=(ANS(1 )-U(N-10)*SIGIKV{4))/2.0 

G(1)=SIGINV(1)*SIGINV<5)+SIG(5)/DET 

G(2)=SIGINV(2)*SIGINV(5)-SIG(3)/<2.0*DET) 

G(3)=SIGINV<3)*SIG1NV<5)-SIG<2)/<2.0*DET) 

G(4)=SIGINV(4)*SIGINV(5) 

G(5)=SIGINV(5)*SIGINV(5)+SIG<1)/<2.0*DET) 

G(6)=SIGIWV<6)»SIGINV(5) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

0T2(N,8)=ANS(1)-U(N-10)»S1GINV(5) 

G<1)=SIGINV(1)»SIGINV<6)SIG(4)/0ET 

G(2)=SIGINV(2)*SIGINV(6)+SIG(2)/DET 

G<3)=SIGINV(3)»SIGINV(6) 

G<4)=SIGINV<4)*SIGINV(6)SIG(1)/DET 

G(5)=SIGINV(5)*SIG1NV(6) 

G(6)=SIGINV(6)*SIGINV<6) 

CALL  MXTMLT(G,MSUM,ANS, 1,6,1) 

OT2<N,9)=(ANS(1)-W(N-10)*S1G1NV(6))/2,0 

RE'^URN 

ELD 


subroutine  DT3VAL(DT3) 

SETS  THE  VALUES  OF  THE  MATRIX  0T3 

INCLUDE  'COM  DEF' 

INCLUDE  'WVEC  DEF' 

REAL  DT3<9,22),U41,U42,U43,U51,U52,U61,U62,U71,U72,U73,ua'!,U82, 
C  U91,U92,U93 
INTEGER  I,J 
DO  100  1=1,3 
DO  200  J=1,22 
DT3(1 ,J)=0.0 
200  CONTINUE 
100  CONTINUE 

DT3(1,15)=-W(8)/(W(5)»»2) 

DT3(1,18)=1.0/W(5) 

DT3(2,15)=W(10)/{W(5)**2) 

DT3(2,20)=1.0/W(5) 

DT3<3,15)=0.5*U(6)/(U(5)**2) 

DT3(3,16)=0.5/W(5) 

U41=R*(2.0*SIGINV(1)*W(9)+SIGINV(2)*W<11)+0.5*SIGINV(3)*(W(1 ) 

C  -W(2)))/(4.0*WC5)) 

U42=W(8)/(2.0*U{5)) 

U43=SIGINV(1 )*MU(1 )*SIGINV(2)»MU(2)+SIGINV(3)*HU(3) 

U4=U41 -U42*U43 
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DT3(4,1)=-SIGINV(1)*UA2 

DT3(4,2)=-SIGINV(2)*U42 

0T3(4,3)=-SIGINV(3)*U42 

U51=R*(SIGINV(1)*U<11)+2.0*SIGINV(2)*U<12)+0.5*SIGINV(3)»(W(3) 

C  -U(4)))/(2.0»W(5)) 

U52=W(10)/U(5) 

U5=U51-U52*U43 

0T3(5,1)=-SIGINV(1)*U52 

0T3(5,2)=-SIGINV(2)*U52 

DT315,3)=-S1GINV(3)»U52 

U61=R*(S!GINV(1)*(W(1)-U(2))+SIGINV(2)»<W(3)-U(4))+SIGINV(3) 

C  *(W(7)*W<5)))/(4.0*U(5)) 

U62=W(6)/(2.0*W(5)) 

U6=U61 •U62»U43 
DT3{6  1)=-SIGINV(1)*U62 
DT3(6,2)=-SIGINV(2)*U62 
DT3(6,3)=-SIGINV(3)*U62 

U71=R»(SIGINV<2)»W(11)+2.0»SIGINV(4)*W<12)+0.5*SIGIN\(5) 

C  •(W(3)-U<4)))/(4.0*U(5)) 

U72=W(10)/(2.0*U(5)) 

U73=SIGINV(2)»MU(1)+SIGIMV(4)*MU(2)+SIGINV(5)»HU(3) 

U7=U71-U72»U73 

DT3(7,1)=-SIGINV(2)*U72 

DT3(7,2)=-SIGINV<4)*U72 

DT3(7,3)=-SIGINV(5)*U72 

U81=R*(SIGINV(2)*(U(1)-W(2))+SIGINV<4)»{U(3)-U(4))+SIGINV(5) 

C  *(W(7)»U(5)))/(4.0*U(5)) 

U82=U(6)/(2.0*U(5)) 

U8=U81 ■U82*U73 
0T3(8,1)=-SIGINV(2)»U82 
DT3<8,2)=-SIGINV(4)»U82 
DT3<8,3)=-SIGINV(5)*U82 

U91=R*{SIGINV(3)*(W(1)-W(2))+SIGINV(5)*(W(3)-W(4))+S!GlNVi'6) 

C  •(W(7)+W(5)))/(8.0»U(5)) 

U92=W(6)/(4.0»W(5)) 

U93*SIGINV<3)*HU(1)+SIGINV(5)*MU(2)+SIGIMV(6)*MU(3) 

U9=U91-U92*U93 

0T3(9,1)=-SIGINV<3)»U92 

0T3(9,2)=-SIGINV(5)*U92 

0T3(9,3)*-SIGINV(6)*U92 

TEMPI =<(R*W(11)/2.0)-MU(2)*W(8))/(2.0*W(5)) 

TEMP2=(0.25*R*(W(1)-W<2))-MU<3)»U(8))/(2.0*U(5)) 

TEMP3  (R*W(9)-MU(1)*U(8))/(2.0*W(5)) 

DT3(4,4)=-SIGINV<1)*U4 

0T3(4,5)=-2.0*SIGINV(2)*U4-(SIG(6)*TEMP1-SIG(5)*T£MP2)/DET 
0T3(4,6)=-2.0*SIGINV(3)*U4(SIG(4)»TEMP2-SIG(5)*TEMP1)/DET 
DT3<4,7)=-SIGINV(4)*U4-(SIG(3)*TEHP2-S!G<6)*TEMP3)/0ET 
0T3(4,8)=-2.0*SIGINV(5)*U4(2.0*SIG(5)»TEMP3-SIG(3)*TEMP1-SIGv?) 
C  *TEMP?)/DET 

0T3{4,9)^  S1G1NV(6)»U4-(SIG(2)»TEMP1-SIG(4)»TEHP3)/0ET 
TEMP1=(R*W(12)-MU(2)*U<10))/W(5) 
TEMP2=(0.25*R*(U(3)U(4))-MU(3)*W(10))/U{5) 
TEMP3=(R*U(11)/2.0-MU(1)*W(10))/U(5) 

0T3(5,4)=-SIGINV(1)*L'5 

DT3(5,5)=-2.0»SIGINV<2)*U5-<SIG(6)*TEMP1SIG(5)*TEMP2)/DET 
0T3(5,6)=-2.0*SIGINV(3)*U5-<SIG(4)»TEMP2-SIG(5)*TEMP1)/DET 
DT3(5,7)=-SIGINV(4)*U5(SIG(3)»TEMP2-SIG{6)»TEMP3)/DET 
DT3(5,8)=-2.0»SIGINV(5)*U5-(2.0»SIG(5)*TEMP3-SIG(3)*T£MP1  $10(2) 
C  •[EMP2)/DET 

DT3(5,9)=SIGINV(6)*U5-(SIG{2)*TEMP1-SIG{4)*TEMP3)/DET 
DT3(7,4)=-SIGINV(1)*U7(SIG(5)»TEMP2SIG(6)*TEMP1)/(2.0*DET1 
OT3<7,5)=-2.0*SIGINV(2)»U7-(SIG(6)*TEMP3SIG{3)*TEHP2)/(2.0‘DET) 
0T3(7,6)=-2.0*SIG1NV(3)*U7-(2.0*SIG(3)*TEMP1-SIG(5)‘TEMh3-S1G(2) 
C  •TEMP2)/(2.0*DET) 

DT3(7,7)  -SIGINV(4)*U7 

DT3(7,8)=-2.0*SIGINV(5)*U7-(SIG(1)»TEMP2SIG(3)*TEMP3)/(2.0*DET) 

DT3(7,9)=SIGINV(6)*U7(SIG(2)*TEMP3SIG(1)»TEMP1)/(2.0»DET) 

TEMP1=(R*(W(3)U(4))/2.0MU(2)*W<6))/(2.0*U(5)) 
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TEMP2=(R*(W(7)+W(5))/2.0-MU(3)*U(6))/(2.0*U(5)) 

TEMP3=<R»<U(1)-W(2))/2.0-MU(1)*W(6))/(2.0*W{5)) 

DT3(6,4)=-SIGINV{1)*U6 

DT3(6,5)=-2.0*SIGINV(2)»U6-(SIG(6)*TEMP1-SIG{5)*TEMP2)/DET 
DT3(6,6)=-2.0»S!GINV(3)*U6-(SIG(4)*TEMP2-SIG(5)*TEHP1)/DET 
DT3(6,7)=-SIGINV{4)*U6-(SIG(3)*TEMP2-SIG{6)*TEHP3)/DET 
OT3(6,8)=-2.0*SIGINV(5)*U6-(2.0*SIG{5)*TEHP3-SIG(3)*TEMP1SIG(2) 
C  *TEHP2)/DET 

DT3(6,9)=-SIGINV(6)*U6(SIG(2)*TEMP1-S1G(4)*TEMP3)/DET 
DT3(8,4)=SIGINV(1)*U8-(SIG(5)*TEMP2-SIG(6)»TEMP1)/0£T 
DT3(8,5)=-2.0*SIGINV(2)»U8-(SIG(6)*TEHP3-SIG(3)*TEMP2>/DET 
DT3(8,6)  =  -2.0*SIGmV(3)»U8-(2.0»S?G(3)*TEMP1-SIG(5)*TEMP3-SIG(2) 
C  •TEMP2)/0ET 

DT3(8,7)=-SIC1NV(4)*U8 

DT3<8,8)=-2.0*SIGINV(5)»U8-(S1G(1)*TEHP2-S1G(3)*TEMP3)/DET 
DT3(8,9)=-S1G1NV(6)*U8-(SIG(2)»TEMP3-S1G(1)*TEMP1)/DET 
0T3(9,4)=-SIGINV{1)*U9-(S1G<5)*TEMP1-SIG(4)*TEMP2)/(2.0»DET) 
DT3(9,5)=-2.0*SIGINV(2)*U9-<2.0*SIG(2)*TEMP2-SIG(3)*TEMP1-SIG(5) 
C  •TEMP3)/(2.0*0ET) 

DT3(9,6)=-2.0*SIG!NV(3)*U9-(SIG(4)*TEMP3-SIG(2)*TEHP1)/(2.0*DET) 

DT3{9,7)=-SIGIhV<4)*U9-<SIG(3)*TEMP3SIG(1)*TEMP2)/(2.0*DET) 

DT3(9,8)=-2.0»SIGINV(5)»U9-(SIG(1)*TEMP1-SIG(2)*fEMP3)/(2.0*DET) 

DT3<9,9)=-SIG1NV(6)*U9 

DT3<4,10)=U41/R 

DT3(5,10)=U51/R 

DT3(6,10)=U61/R 

DT3(7,10)=U71/R 

DT3(8,10)=U81/R 

OT3(9,10)=U91/R 

00  300  1=4,9 

00  400  J=11,22 
OT3<I,J)=0.0 
400  CONTINUE 
300  continue 

OT3<4,11)=<R*SIGINV(3))/<8.0*W<5)) 

DT3<4,12)=-0T3(4,11) 

0T3(4,15)=-U4/U(5) 

OT3<4,18)=-U43/(2.0*W(5)) 

0T3<4,19)=(R*SIGINV<1))/<2.0*U(5)) 

OT3<4,21)=(R*SIGINV<2))/(4.0*W<5)) 

OT3<5,13)=(R*SIGINV<3))/(4.0»U(5)) 

0T3(5,14)=DT3(5,13) 

OT3{5,15)=U5/W(5) 

OT3(5,20)=-U43/W(5) 

0T3<5,21)=0T3(4,19) 

0T3(5,22)=(R»SIGINV{2))/U<5) 

0T3(6,11)=(R*SIGINV(1))/(4.0*W<5)) 

0T3(6,12)=DT3(6,11) 

DT3(6,13)=0T3(4,21) 

0T3(6,14)=0T3(6,13) 

0T3(6,17)=(R*SIGINV<3))/(4.0*W(5)) 

0T3<6,15)=-U6/W(5)+0T3(6,17) 

0T3(6,16)=-U43/(2,0*W(5)) 

0T3(7,13)=(R«SIGINV(5))/<8.0*W(5)) 

0T3(7,14)=-0T3(7,13) 

0T3(7,15)=-U7/W(5) 

DT3(7,20)=-U73/(2.0*W(5)) 

0T3(7,21)=DT3(4,21) 

DT3(7,22)=<R*SIGINV(4))/(2.0*W(5)) 

0T3(8,11)=DT3(4,21) 

0^3(3, :2)=  0T3(4,21) 

DT3<8,13)=0T3(7,22)/2.0 

DT3(8,H)=-DT3<8,13) 

OT3(8,15)=-U8/W<5)+OT3(7,13)*2.0 

0T3(8,16)=0T3(7,20) 

DT3(8,17)=DT3(7,13)*2.0 

DT3(9,11)=DT3(4,11) 

DT3(9,12)=-DT3(9,11) 
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0T3(9,13)=0T3(7,13) 

0T3(9,H)=-DT3{9,13) 

DT3{9,17)=(R*SIGINV(6))/(8.0*U(5)) 

DT3{9,15)=-U9/W(5)+DT3(9,17) 

DT3(9,16)=U93/(4.0»U(5)) 

RETURN 

END 


SUBROUTINE  MAKES 

*  CREATES  THE  MATRICES  S  AND  ST 

INCLUDE  'PMUCOM  DEE' 

INTEGER  I.J 
DO  10  1=1,9 
DO  20  J=1,6 
SCI,J)=0.0 
20  CONTINUE 
10  CONTINUE 
S(1,1)=1.0 
S(2,2)=0.5 
S(3,3)=0.5 
S(A,2)=0.5 
S(5,4)=1.0 
S(6,5)=0.5 
S(7,3)=0.5 
S(8,5)=0.5 
S(9,6)=1.0 
DO  30  1=1,6 
DO  40  J=1,9 

ST(I,J)=S(J,I) 

40  CONTINUE 
30  CONTINUE 
RETURN 
END 


SUBROUTINE  PMUS( INP,OUTP1 ,OUTP2) 

•  COMPUTES  THE  VALUES  OF  PMU  AND  PSIG 

INCLUDE  'PMUCOM  DEF' 

INCLUDE  'COM  DEF' 

REAL  INP(3,3),0UTP1(3,3),PR00(9,9),TEMP(9,6),0UTP2(6,6) 
INTEGER  I,J 

CALL  TENSR<INP,INP,PR00,3,3,3,3) 

CALL  MXTMLT(PRCX),S,TEMP,9,9,6) 

CALL  MXTMLT(ST,TEMP,CIUTP2,6,9,6) 

DO  10  1=1,6 
DO  20  J=1,6 

0UTP2(I , J)=(2.0*OUTP2(I , J))/(REAL(NUMTR)*NUMSIM) 
20  CONTINUE 

10  CONTINUE 
DO  30  1=1,3 
DO  40  J=1,3 

OUTP1(I,J)=INP<I,J)/REAL(NUMTR) 

40  CONTINUE 

30  CONTINUE 
RETURN 
END 
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SUBROUTINE  TENSR(A,B,C,A1 ,A2,B1 ,B2) 

»  COMPUTES  THE  TENSOR  PRODUCT  OF  A  AND  B  OUPTUT  TO  C 

INTEGER  I,J,IC,L,A1,A2,B1,B2 

REAL  A(A1 ,A2),B(B1 ,B2),C(A1*B1 ,A2*B2) 

DO  10  1=1, A1 
DO  20  J=1,A2 
DO  30  K=1,B1 
00  40  L=1,B2 

C(K+(I-1)*B1,L+(J-1)*B2)=A(!,J)*B(K,L) 
40  CONTINUE 

30  CONTINUE 

20  CONTINUE 
10  CONTINUE 
RETURN 
END 


SUBROUT  I NE  BNOEST ( SEPM , PMU , PS I G , SEPMU , SEPS I G , C I ) 

*  CALCULATES  THE  FIRST  AND  SECOND  ORDER  VARIANCE  ESTIMATES  AND 

*  CORRESPONDING  CONFIDENCE  INTERVAL  SIZES 

INTEGER  I.J 

REAL  TEMP1(3,3),TEMP2(6,6),TEMP3(3,6),TEMP4(6,3),SSQMU<9), 

C  SS0SIG(36),SSQMSG<18),SSQSGM(18),SEPM{9,9),PSIGSQ(36,36), 

C  PSIG(6,6),SEPMU<3),SEPSIG(6),TEMP(36),PMUSQ(9,9),PMU(3,3), 

C  PMUSIG(18,18),PSIGMU(18,18),CI(2).STPHU(9).STPSIG(36) 

C  ,BND,BN01,8ND2,BND3 

DO  10  1=1,3 
DO  20  J=1,3 

TEMP1(I,J)=SEPM(I,J) 

20  CONTINUE 

DO  30  J=1,6 

TEMP3<I,J)=sePM(I,J+3) 

30  CONTINUE 

10  CONTINUE 

CALL  STRING(TEMP1,SSOMU,3,3) 

CALL  STRING<TEMP3,SSQMSG,3,6) 

DO  40  1=1,6 
DO  50  J=1,3 

TEMP4(I,J)=SEPM(I+3,J) 

50  CONTINUE 

DO  60  J=1,6 

TEMP2{I ,J)=SEPM<I+3,J+3) 

60  CONTINUE 

40  CONTINUE 

CALL  STRINGCTEMP4,SSOSGM,6,3) 

CALL  STRING(TEMP2,SSQSIG,6,6) 

CALL  MXTMLT(PMU,SEPMU,TEMP,3,3,1) 

CAL L  MXTML T ( SEPMU , TEMP , BND , 1 , 3 , 1 ) 

CALL  MXTMLT(PSIG,SEPSIG,TEMP,6,6,1) 

CALL  MXTMLT(SEPSIG,TEMP,BN01, 1,6,1) 

BND1=8ND1+BN0 

CALL  TENSR(PMU,PMU,PMUSQ,3,3,3,3) 

CALL  MXTML T ( PMUSO , SSOMU , TEMP , 9 , 9 , 1 ) 

CALL  MXTMLTCSSOMU, TEMP, BND, 1 ,9,1 ) 

CALL  TENSR(PSIG,PSIG,PSIGS0,6,6,6,6) 

CALL  MXTMLT(PSIGS0,SSQSIG,TEMP,36,36,1) 

CALL  MXTMLTCSSOSIG, TEMP, 8ND2, 1,36,1) 

BND2=(BN0+BN02)/'2.0 

CALL  TENSR(PMU,PSIG,PMUSIG,3,3,6,6) 

CALL  MXTMLT(PMUSIG,SS0MSG, TEMP, 18, 18,1) 

CALL  MXTML T ( SSOMSG , TEMP , BND , 1 , 18 , 1 ) 

CALL  TENSR(PSIG,PMU,PSIGMU,6,6,3,3) 

CALL  MXTMLT(PSIGMU,SSQSGM, TEMP, 18, 18,1) 

CALL  MXTMLT(SSQSGM, TEMP, BND3, 1,18,1) 
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BND2=BND2'»(  (  BND+BND3  )/2 . 0 ) 

CALL  STRING(PMU,STPMU,3,3) 

CALL  STRING(PSIG,STPS1G,6,6> 

CALL  MXTMLT(SSQMU,STPMU,BND, 1,9,1) 

CALL  MXTMLT ( SSQS I G , STPS I G , BND3 ,1,36,1) 

BN02=BND2+(((BNO+BND3)»*2)/A.O) 

BN0=BND1+BND2 

CI(1)=1.96*SQRT(BND1) 

CI<2)=1.96*SQRT(BND) 

RETURN 

END 


SUBROUTINE  STRINGCAI , A2,N1 ,N2) 

*  STRINGS  CUT  THE  MATRIX  A1  AS  A  VECTOR  A2 

REAL  A1(N1,N2),A2(N1*N2) 

INTEGER  I,N1,N2,J 
DO  10  1=1, N1 
DO  20  J=1,N2 

A2((1-1)*N2+J)=A1(1,J) 

20  CONTINUE 
10  CONTINUE 
RETURN 
END 


REAL  FUNCTION  HX(X,Y) 

REAL  X,Y,SPHER(3),TEMP 
INCLUDE  'COM  DEF' 

SPHER(1)=R*SIN(X)*COS(Y)-MU(1) 

SPHER(2)=R*SIN(X)»SIN{Y)-MU(2) 

SPHER(3)=R*COS<X)MU<3) 

TEMP=-0.5*<(<SPHER<1)*»2)*SIGINV(1))+(2.0*SPHER(1)*SPHER(2) 

C  *SIGINV(2))»<2.0*SPHER(1)*SPHER(3)»SIGINV(3)) 

C  ♦((SPHER(2)**2)*SIGINV(A))+(2.0*SPHER{2)»SPHER(3) 

C  *SIGINV<5))+<<SPHER(3)**2)*SIGINV(6))) 

IF  (TEMP  .GT.  -27.0  )  THEN 
HX=EXP(TEMP)/DENOM 
ELSE 

HX=0.0 
END  IF 
RETURN 
END 

REAL  FUNCTION  CC1111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC1111=COS(X)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1311(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC1311=COS(3.0*X)«COS(Y)»HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CC3111(X,Y) 

EXTERNAL  HX 
REAL  X,v 

CC3111=(C0S(X)**3)*C0S(T)*HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CS1111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1111=C0S(X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1311(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1311=COS(3.0«X)*SIN{Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS3111={COS<X)»*3)*SIH(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1100(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1100=SIN<X)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1200<X,Y) 
EXTERNAL  HX 
REAL  X,Y 

SC1200=SIN(2.0*X)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1300<X,Y) 
EXTERNAL  HX 
REAL  X,Y 

SC1300=SIN(3.0*X)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  SC2111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC2111=(SIN(X)*»2)*COS<Y)»HX(X,Y) 

RETURN 

END 

real  function  SC3121(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC3121=(SIN(X)**3)*(C0S(Y)»*2)»HX<X,Y) 

RETURN 

END 

REAL  FUNCTION  SS2111(X,Y) 
external  HX 
REAL  X,Y 

SS2111=(SIN(X)**2)«SIH<Y)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  SS3112(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS3112=(SIN(X)**3)*S1N(2.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS3121(X,Y) 

EXTERNAL  HX 
REAL  X.Y 

SS3121=(SIN(X)»*3)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CC1113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC1113=COS(X)»COS(3.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1131(X,Y) 

EXTERNAL  HX 
REAL  X.Y 

CC1131=C0S(X)*(C0S<Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1211(X,Y) 

EXTERNAL  HX 
REAL  X.Y 

CC1211=COS(2.0*X)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1331(X,Y) 

EXTERNAL  HX 
REAL  X.Y 

CC1331*COS(3.0*X)»<COS(Y)»*3)*HX(X,Y) 

RETURN 

END 

HEAL  FUNCTION  CCU11(X,Y) 

EXTERNAL  HX 
REAL  X.* 

CCU11  )S(4.0*X)»COS<Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC1511(X,Y) 

EXTERNAL  HX 
REAL  X.Y 

CC1511=COS<5.0*X)*COS(Y)*HX(X,T) 

RETURN 

END 

REAL  FUNCTION  CC1531(X,Y) 
external  HX 
REAL  X.Y 

CC1531=COS(5.0*X)*(COS(Y)*'‘5)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CC2111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC2111=(COS(X)»*2)*COS(Y)*HX{X,Y) 

RETURN 

END 

REAL  FUNCTION  CC3113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC3113=(COS(X)»*3)*COS(3.0*Y)»HX<X,Y) 

RETURN 

END 

REAL  FUNCTION  CC3131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC3131=<C0S{X)»»3)*{C0S(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC5111<X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC5111=(COS<X)**5)*COS(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CC5113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC5113*(COS(X)»*5)*COS(3.0*Y)»HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CC5131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CC5131=(COS(X)»*5)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CS1113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1113=C0S(X)*SIN(3.0»Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1131=COS<X)*(SIN(Y)*«3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1211(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1211=COS(2.0*X)*SIN(Y)»HX(X,Y) 

RETURN 

END 
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REAL  FUNCTION  CS1331(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1331=COS(3.0*X)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CSU11(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CSK11=COS(4.0»X)»SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1511(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1511=COS(5.0*X)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS1531vX,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS1531=COS(5.0*X)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS2111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS2111=(COS(X)**2)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3113<X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS3113*<COS<X)**3)*SIN(3.0»Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS3131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS3131=(COS(X)**3)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  CS5in(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS5111=(C0S(X)»»5)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTIO-I  CS5113(X,Y) 
external  HX 
REAL  X,T 

CS5l’3=(COS(X)»*5)»SIN(3.0»Y)»HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  CS5131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

CS5131=<COS(X)»»5)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  SC1121(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1121=SINvX)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1221(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1221=SIN{2.0*X)*(COS<Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1321(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1321=SIH(3.0*X)»<COS<Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1400(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SCH00=SIN(4.0*X)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1421(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1421=SIN(4.0*X)*<COS(Y)**2)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1500(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1500=SIN(5.0*X)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1521<X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1521=SIN(5.0*X)*(COS(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC3100(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC3100=(SIH(X)*«3)*HX(X,Y) 

RETURN 

END 


REAL  FUNCTIC3N  SC4111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC4111=(SIN(X)**4)*C0S(Y)*HX{X,Y) 

RETURN 

END 

REAL  FUNCTION  SC4113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC4113=(SIN(X)»»4)*COS(3.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC4131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC4131=(SIN(X)**4)*(COS(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5112(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC5112=(SIN(X)**5)*COS(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5114(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC51U=(SIN(X)**5)*COS<4.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC512UX,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC5121=(S1N<X)*»5)*(C0S(Y)**2)*HX<X,Y) 

RETURN 

END 

REAL  FUNCTION  SC5141{X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC5141=(SIH(X)**5)*(COS(Y)**4)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1212(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SC1212=SIN(2.0*X)«COS(2.0»Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SC1412(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SCU12=SIN(4.0*X)*COS(2.0*Y)*HX<X,Y) 

RETURN 

END 
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REAL  FUNCTION  SS1112(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1112=SIN(X)»SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1121(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1121=SIN(X)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1212(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1212=SIN(2.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1221(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1221=SIN(2.0*X)*(SIN<Y)»*2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1312(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1312=SIN<3.0*X)»SIN(2.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1321<X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1321=SIN(3.0*X)*<SIN(Y>»*2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SSU12(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1412=SIN(4.0*X)*SIN(2.0*Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SSU21{X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SSU21=SIN(4.0*X)»(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS1512{X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1512=SIN{5.0*X)*SIN(2.0*Y)»HX(X,Y) 

RETURN 

END 


REAL  FUNCTION  SS1521(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS1521=SIN(j.0*X)»(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS4111(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS4111=(SIN(X)**4)*SIN(Y)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS4113(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS4113=<SIN(X)**4)*SIN(3.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS4131(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS4131=(SIN(X)»*4)*(SIN(Y)**3)*HX(X,Y) 

RETURN 

END 

REAL  FUNC: ION  SS5112rX,Y) 

EXTERNAL  HX 
REAL  X.Y 

SS5112=(SIN(X)**5)»SIN(2.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS5114(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS5114s(SIN(X)**5)*SIN(4.0*Y)»HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS5121(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS5121=(SIN(X)**5)*(SIN(Y)**2)*HX(X,Y) 

RETURN 

END 

REAL  FUNCTION  SS3141(X,Y) 

EXTERNAL  HX 
REAL  X,Y 

SS5141=(SIN(X)**5)*(SIN(Y)**4)*HX(X,Y) 

RETURN 

END 


t 


COM  DEF 

REAL  R,SIG(6),SIGINV(6),DET,DENOH,MU(3),PI,A,B,ERRABS,ERRREL 
C  ,MSIG(3,3),MUSET(3),SIGSET(3,3),SEPSET 
INTEGER  IRULE,ISEED,NLMTR 

COMMON/COM1/R,SIG,SIGINV,DET,OENOM,MU,ERRABS,ERRREL,MSIG, 

C  MUSET , SI GSET , SEPSET , I  RULE , NUNTR , I  SEED 

PARAMETER(PI=3.  U1S92(') 
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WVEC  DEF 

REAL  U(12),CC{17),CS(18),SC(19),SS(18) 
COMMON/UCC/W, CC, CS, SC,  SS 


PMUCOM  DEF 

REAL  S(9,6),ST(6,9),NUMSIM 
COHMON/COM3/S , ST , NUMS I M 


SEPCC3M  DEF 

REAL  CXX , CYY , CZZ , SQP I! NV, QUP I TU, CSAN T , C , RAD  1 US ,MN(3 ) 
COMMON/COM2/CXX , CYY , CZZ , CSANT . C , RAD  1 US ,MN 
PARAHE  T  ER ( SQP 1 1 N V=0 . 564 1 896, QUP I T W=0 . 05626977) 
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